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COMPUTER  ROUTINES  FOR  USE  IN  SECOND-ORDER 
V01.TERRA  MODELING  OF  EMI 


I.  INTRODUCTION 

Research  described  here  consists  of  the  development  of  certain  computer 
routines  to  aid  Volterra  modeling  of  weakly  nonlinear  systems  [1],  (2).  Volterra 
series  representation,  a  dynamic  generalization  of  the  familiar  power  series,  is 
ideal  for  representing  devices  and  systems  with  frequency-dependent  mild  non- 
linearities  as  in  the  case  of  a  transistor.  The  technique  has  been  applied  to 
the  analysis  of  communication  receiver  response  to  radio  frequency  interference[31 . 
Other  applications  include  Intermodulation  distortion  analysis  of  transistor 
feedback  amplifiers  (4),  nonlinear  characterization  of  IMPATT  diodes  and  micro- 
wave  amplifiers  (51,  and  analysis  of  channels  with  soft  limiter. 

Rome  Air  Dev«.lopment  Center  has  in  the  recent  past  supported  several 
efforts  to  put  this  analytical  tool  to  use  in  the  electromagnetic  interference 
and  compatibility  field.  (In  practical  terms  one  of  tl»e  major  outcomes  of  the 
efforts  is  the  lAP  program,  a  computer  program  for  the  prediction  of  Intra- 
System  Electromagnetic  Compatibility).  A  current  direction  of  interest  is  the 
estimation  of  Volterra  kernels  of  a  system  from  its  experimentally  observed 
input-output  responses.  Interest  in  this  black-box  approach  arises  for  several 
reasons,  the  most  salient  being  cost  effectiveness  in  testing,  and  simplicity 
of  resulting  models  for  complex  on-board  communication  systems. 

Weiner  and  Ewen  (I),  (2]  have  provided  an  approach  to  finding  the  para¬ 
meters  of  the  kernels,  specifically  the  poles  and  residues  of  the  multivariable 
transfer  functions  H  (sp....,s^).  The  poles  and  residues  of  the  linear  TF, 

H.(s),  are  determined  using  Jain's  method  [6],  [ 7) (pencil-of-functlons  identi¬ 
fication  method).  Then,  for  somewhat  larger  amplitudes,  where  the  quadratic 
TF  H2(s,,S2)  has  non-negligible  influence,  the  contribution  y,(t)  is  determined 
by  subtracting  y^Ct),  the  predicted  response  of  Hj(s),  from  yft).  The  poles 
of  the  quadratic  TF  are  known  in  terms  of  the  poles  of  the  linear  TF,  so  that 
only  the  residues  of  H,(Sj,  s.,)  need  be  —  and  in  principle,  can  be  —  deter¬ 
mined.  A  similar  procedure  is  adopted  for  determining  the  parameters  of  the 
cubic  kernel,  and  so  on. 

The  computer  programs  presented  in  this  report  are; 

IGRAM  Program  for  black-box  identification  of  a  linear  transfer 
function  using  pencil-of-functlons  method. 

HPMINV  High  precision  matrix  inversion  routine  for  use  in  deter¬ 
mination  of  the  residues  of  the  quadratic  Volterra  TF. 

Perturbation  theory  and  iterative  corrections  are  used  to 
enable  accurate  Inversion  even  for  wideband -system  matrices. 


1 


II.  PENCll.  OF  FUNCT10N.‘:  METHOD  FOR  IDENTIFIC.^TION 
OF  NETWORK  TR^VN.RFER  FUNCTIONS 


Determining  the  moilcl  of  .t  network  from  its  observed  input-output  responses 
represents  the  inverse  of  the  analysis  problem.  Interest  in  this  arises  from  the 
frequent  need  for  a  relatively  simple  mathematical  description  of  the  system  so 
that  behavior  for  otlier  anticipated  inputs  may  be  predicted  up  to  acceptable 
accuracies.  Like  the  analysis  problem,  there  are  several  approaches  available 
in  the  literature  for  the  inverse,  or,  as  it  is  often  called,  the  "identification" 
problem  .  To  name  a  few,  (a)  Prony's  method  (Si,  (b)  gradient  methods,  such 
as  Newton  (9I  and  qua  '  »-l inearization  110j,(c)  least-squares  and  generalized 
least-squares  methods  1 1 1 ) , 1 1 2 ) , (d)  maximum- likelihood  methods  1 13) , [ 14] , [ 15] . 

All  of  the  methods  stated  above  possess  certain  advantages  and,  as  may  be 
expected,  certain  disadvantages  peculiar  to  each  partl''ular  method.  Stated 
very  broadly,  sensitivity  to  noise,  slow  convergence  to  the  solution,  and  ex¬ 
cessive  computational  complexity  are  some  of  the  possible  disadvantages.  The 
objective  of  this  section  is  to  describe  in  a  seml-rlgorous  way  the  pencil-of- 
functions  identification  method  [7].  Further,  in  this  section  the  method  is 
extended  to  the  case  where  general  first  ordei  -s,  U^fs)  =  (b^s  +  Cj)/(s+a^) 

are  used  in  the  processing  system  instead  of  ide.  integrators  (note  that  the 
ideal  integrator  is  a  special  case  of  this  filter;  set  b^  =  a^  =  0).  .A  high 
accuracy  FORTRAN  program,  developed  for  the  case  of  ideal-integrator  processing 
units,  is  also  presented. 

The  method  offers  the  advantages  of  mathematical  simplicity,  closed-form 
solution  to  the  problem,  which  is  optimal  in  the  generalized  least-squares  sense 
and  suboptiraal  in  the  strict  least-squares  sense,  and  robustness  to  noise.  The 
disadvantage  of  the  method  is  that  the  variances  of  additive  noise,  when  present, 
must  be  determined  in  a  separate  experiment  in  order  that  unbiased  parameter 
estimates  may  be  computed. 


2.1  THEORY 


The  problem  of  identifying  the  transfer  function  of  a  network  from  its  innut- 
output  responses  can  be  formulated  in  discrete-time  domain  as  follows.  Given  the 
pair  x(k),  y(k)  (or  noise-corrupted  versions  of  these;  call  them  Tj(k),  v(k))  find 
the  transfer  function  K(2)  that  produces  a  response  matching  y(k)  (or  v(k))  when 
excited  by  x(k)  (or  u(k)).  More  specifically,  this  involves  determination  of  the 
parameters  a^^,  b^  in 

Y(z)  =  H(2)  X(2) 

b  ■^'bz  ■f‘.**"l‘bz 

=  (-2 - - s—  )  x(2)  (1) 

1  +  a, z  +...+  a  z 
1  n 

from  experimental  input-output  data.  Note  (1)  can  be  written  in  time  domain  as 
n  n 

y(k)  =  -  E  a  y(k-i)  +  E  b^  x(k-i),  y(k)  =  0  for  k  <  0  (2) 

i=l  i«0  ■ 

A.  Measurement  Signals 

Before  proceeding  with  the  solution  of  the  problem  via  pencil-of-function 
method,  let  us  note  that  (2)  can  be  written  as 

[y(k)  y(k-l)  .  .  .  y(k-n)  -x(k)  -x(k-l)  .  .  .  -x(k-n)]  0=0  (3) 

or,  more  concisely, 

f  (k)"^  0  =  0  (4) 

where  the  2n  +  2  dimensional  vector  £(k)  has  the  obvious  definition  and  the  vector 
also  2n  +  2  dimensional,  is  given  as 

(5) 

Equation  (4)  represents  a  geometrical  constraint  upon  the  vectors  ^(k) ,  namely 
that  they  are  all  orthogonal  to 

To  cast  the  problem  of  identifying  ^  into  a  generalized  least-squares 
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problem  consider  the  measurement  system  of  Fig.  1.  The  matter  of  choosing  the 


first  order  filters  which  make  up  the  measurement  units  M^(2)  will  not  be 
considered  here.  It  will  suffice  to  say  they  must  be  chosen  so  that  the 
cuti  ff  frequencies  of  Mj<r)  span  the  frequency  range  of  interest  (l.e.  they  are 
spread  in  the  pass-band  of  the  network  under  test).  Kegardless  of  the 
choice  of  the  measurement  units  the  following  useful  observation  arises. 

Let 

0{k)  *  {y  (k)  y  (k)  ...  y  (k)  -x  (k)  -x, (k)  ...  -x  (k))^ 

—  oi  n  o  I  n 

and  deline  the  matrix 


*^00 

*^01  ■ 

i  C 

!  10 

f 

*^11 

‘^In 

1 

t  “^nO 

*^01 

c 

nn 

where  c^^  are  the  coefficients  of  the  polynomial 


n 

I 

i-0 


ji 


-1  ^  -1  ”  -1 
-  >■  '  (1-P.z  ‘)  -  (I-Q.2  *) 

i=l  i=i+l 


(k) 

y„(k) 

yj(k)  ii- 

1 

N 

1 

.',(k) 

1-P 

n 

V  (k) 

"  n 

*!  -1 1 
ll-QlS  1 

t- 

L. 

-Q,2  : 

-i» 

;Ck) 

i-Pjz-i; 

-  ^  -  —  —  \ 

Xj^(k) 

l-P^z”^  ■ 

s^(k) 

1-P  2*^ 
n 

X  (k) 

n 

[l-Ql-'*' 

1 

l-Q  z"^ 
L  n  _ 

^ - ^  - 

- 

J 

M,{2) 


MjU) 

Measurement  vector  Jy^Ck) - y^(k)  -x^Ck)  ....  -Sj^(k)l’‘  =  5(k) 

First  order  filter.s  U.Cz)  =  (l-P.z-h/d-Q.z'S 


Fig.  1.  Measurement  svstem 


where 


On  Che  other  hand  the  measurement  signals  have  the  following  represen¬ 


tation.  Consider  the  output  measurement  signal  y^(k);  its  transform  is 


Y,(e)  =  M,(e)Y(z) 


i  1-P„z 

=  n  — - 


-1 


i=l  l-Qj^z 


zi 


D( 


~  n  (i-p  z~^)  n  (i-Q  z"^)Y(z) 

£=1  ^  £=i+l 


D(z)  ^‘^oi  ‘^li 


Cnil  i  Y(z) 


where 


D(z) 


=  n  (1  -  Q  z  ") 
£=1  ^ 


and  the  numbers  are  the  same  as  defined  in  equation  (6).  Now  the  output 
measurement  vectors  may  be  written  as 


Y(z)  =  (Yq(2)  . 


where  the  (n+1)  x  (n+1)  matrix  C  =•  given  by  the  numbers  defined 

above.  In  a  like  manner,  we  have 

X(z)  =  [Xq(z)  .  .  •X„(z)]'^  =  ^j  c'^CxCz) 

We  can  now  state  and  prove  the  measurement  filter  theorem. 


Theorem  -  If  the  signals  y(k)  and  x(k)  satisfy  (8)  for  some  parameter 
vector  (a,  M,  then  the  measurement  vectors  satisfy  the  orthogonality 
condition 

=  0  for  all  k 


[a  b'^] 


iL(k) 

-x(k) 


whore 


a  =  C”^a 


C“^b 


(9) 


(10) 


(11) 


(12) 


(13) 


(14) 

(15) 


Proof .  The  matrix  C  can  be  shown  to  be  nonsingular.  Therefore  we  can  rewrite 


;>r. 


5  .  . 


(C  "a)*  Y(8) 


-  (C  "b)  X(r.) 


or 


s'"  3(3  m)  " 


Upon  svibscltutlng  (11)  and  (12)  this  oquatlon  yi«'lds 

«  0 


1(3) 

-x(«) 


The  rusv»lt  sought  by  the  theorem  is  obtaiiuul  immediate.lv  upou  taking  the  inverse 


transform. 
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2.2  APPLICATION  EXAMPLES 


EXAMPLE  1 


As  a  first  example  of  the  Identification  routine,  data  from  the  transfer 
function 

H  (s)  -  -iio!!-.  +  + _ LLo!i _ + _ Opi _ 

s+7(10^)  s+7a0^)  a+(10+J25)10^  s+(lO-J25)  10^ 

was  obtained.  The  input  driving  function  was  a  square-wave  followed  by  a  de¬ 
caying  exponential.  ITjo  cycles  of  square  wave  with  period  0.05hsec  were  used, 
followed  by  the  decaying  exponential  with  time  constant  approximately  equal  to 
O.OSpsec.  This  input  was  selected  based  on  apriori  knowledge  of  the  network 
behavior.  The  spectral  content  of  this  input  should  supply  a  sufficient  amount 
of  energy  to  the  fast  mode  (s  «  -7  x  10^)  and  to  the  slow  and  oscillatory  modes 
for  accurate  identification. 

Five  hundred  points  of  data  were  used  (MP1=500)  with  a  sampling  Interval 
A  “  0.002visec.  The  option  IREM  =1  was  used  in  this  example  because  direct  trans¬ 
mission  could  not  be  assumed.  A  fourth  order  model* (N=A)  was  desired  for  this 
networ*'.  The  identified  poles  and  residues  of  the  model  are  given  below,  to¬ 
gether  with  the  actual  (H^(s))  poles  and  residues. 


H, (s)  Poles 

-7.xl0^’ 

-7.xl0^ 


Identified  Poles  H, (s)  Residues  Identified  Residues 
A 


-7.0x10 

-7.0xlo’ 


-(10.+J25)10  -(10.0+j25.0)l0' 

-(10.-j25)10‘’  -(10.0-j25.0)10^ 


1.0(10  ) 
30.0(10^) 
1.0(10^’) 
1.0(10^’) 


1.0(10”) 

30.0(10^) 

i.o(k/’) 

l.OdO*’) 


As  seen,  the  identification  of  the  transfer  function  was  very  accurate,  and 
the  corresponding  mean-square  and  root-mean  square  errors  for  this  Identification 
are  both  0.0%.  Plots  of  tlie  input-output  data  and  tlie  actual  and  model  responses 
are  shown  in  Fig.  2. 

*  A  rational  transfer  function  with  an  nth  degree  denominator  is  referred  hero  as 
being  of  nth  order. 
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Example  2 

We  examine  the  applicability  of  the  identification  technique  to  responses 
obtained  from  a  wide-band  transistor  amplifier  circuit.  The  schematic  and 
equivalent  model  of  the  circuit  are  shown  in  Fig.  3. 


(b) 

Fig.  3.  (a)  Schematic  of  common  emitter  amplifier  circuit 

(b)  Equivalent  circuit  model 

As  shown  in  Appendix  C,  the  network  transfer  function  is 


H(s)  = 


V2(s) 


_ 8(10^)  s^s-8000  (10*^^)) _ 

(s+, 033 ( 10^) )  (s+.  080 { 10^) )  (s+25 . 2  ( lO*’) )  (s+l205 .  1  ( 10^) ) 


(17) 


The  network  function  can  be  identified  successfully  only  by  performing 
separate  tests  in  tliree  different  frequency  regions: 

(i)  Low  frequency  region  (L) 

(li)  Mid  to  High  frequency  transition  (MH) 

(ill)  Higli  Frequency  region  (H) 

A  discussion  of  these  three  regions  is  given  in  Appendix  C.  Here  we  shall  focus 
primarily  on  the  results  of  identification. 

(i)  Low  frequency  region 

An  adequate  low  frequency  description  of  eq.  (17)  is  given  by  * 

2 

(18) 


Hjfs) 


_ (21.04)  s 

(s+,  033  (i0§’)(s+.  080(10^)) 

It  is  thcref ore  desirable  that  we  seek  .a  second-order  (N=2)  model 

of  the  network  given  by  eq.  (17).  The  input  used  is  a  single  triangular  pulse 
of  duration  125iisec.  One  thousand  points  (Ml’l=1000)  of  input-output  responses, 


*  In  practical  applications,  such  approximations  will  of  course  not  bo  available. 
But  the  designer,  or  even  the  test  engineer,  frequently  has  some  idea  ol  the 
critical  frequencies  of  the  system. 
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sampled  at  A  =  0.25ys  are  used  for  modeling.  The  option  IREM  =  0  was  used 
because  our  low-frequency  model  will  exhibit  direct  transmission. 

The  computer  program  IGRAM  yields  the  following  low-frequency  model: 


H  (s)  =  -  20.125  (s-0. 0015(10^)) (s-K). 0012 (10^)) 
^  (s+0. 034(10^)) (s+0. 075(10^)) 


(19) 


Comparison  of  the  identified  model,  eq.  (19),  with  our  low  frequency  approximation 
(eq.  (18))  shows  close  agreement.  The  rms  error  between  the  measured  network  time 
response  and  the  model  response  is  1.206%.  Plots  of  the  input-output  signals  and 
of  the  actual  and  model  responses  are  given  in  Fig.  4. 

11)  Mid  to  High  frequency  transition 

As  discussed  in  Appendix  C,  an  adequate  mid  to  high  frequency  transition 
description  of  eq.  (17)  is  given  by 


u  (3)=  - 

™  (s  +  25.2(10'’)) 


(20) 


which  is  a  single  pole  function.  Thus,  we  will  attempt  to  model  the  circuit 
with  a  first-order  transfer  function  (N=l).  Since  this  approximate  description 
(eq.  (20))  does  not  exhibit  direct  transmission,  IREM  4  0  should  be  used  in  the 
program.  For  this  Identification,  IREM  =  1  was  selected.  For  reasons  mentioned 
in  Appendix  C,  bias  was  assumed  present  (IBIAS  =  1)  on  the  data.  Five  Iiundred 
points  (MPl  =  500)  of  input-output  signals,  sampled  at  A  =  O.Olpsec,  were  used. 
The  excitation  used  was  a  narrow  band  signal  with  a  center  frequency 
near  the  critical  frequency  (s=-25.2(10  )). 

The  computer  program  IGRAM  yields  the  following  mid  to  high  frequency 
transition  model: 


n  -  _520.-.0J.LQ 

(s+24. 92(10  )) 


(21) 


The  rms  error  between  tills  model  a.id  the  actual  i  etwork  respon  (obtained  from 
the  original  4th  order  transfer  function)  is  1.99;,%.  A  graphical  comparison  is 
given  in  Fig.  5. 
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iii)  High  frequency  region 

An  approximate  high  frequency  description  of  eq.  (17)  is  shown  (Appendix  C) 

to  be 


8(10^)  (s-8000(10^)) _ 

( s+2  5 . 2 ( 10^ ) ) ( s+1205 .1(10^)) 


which  is  a  two-pole  function.  We  will  therefore  attempt  to  model  the 

circuit  with  a  second-order  transfer  function  (N“2).  Direct  transmission 
cannot  be  assumed,  so  that  IREM  =  1  was  used  for  the  model.  Five  hundred 
points  (MPl  =  500)  of  input-output  signals,  sampled  at  A  =  0.00025lisec  were 
used.  A  narrow-band  signal  with  a  center  frequency  of  6300  Mrad/sec  was 
used  for  network  excitation. 

The  computer  program  yields  the  following  high  frequency  model: 


2.79(10  )  ((8-19060(10°)) 
(s+25. 7(10^)) (8+1139.5(10^)) 


The  rms  error  between  this  model  and  the  actvial  network  response  is  0.919%. 

A  graphical  comparison  is  given  in  Fig.  6. 

We  have  now  obtained  a  description  of  the  networ  behavior  in  the  three 
frequency  regions  of  interest.  The  models  obtained  may  be  pieced  together  in  such 
a  way  as  to  synthosize  the  overall  behavior  of  the  .system.  Omitting  the  details 
Involved,  the  following  model  m.ay  be  obtained  from  eq's.  (19),  (21)  and  (23): 


H(s)  = 


62.2(10^)  (s- . 0015(10^) ) (S+.0012 ( 10^) ) (s- 199060 ( 10^) ) 
(s+.034(10^))(s+.075(10^))(.s+2  4.9  (10*-'))  (s+11 39. 5(10^') ) 


Comparison  of  this  model  with  cq.  (17)  shows  favorable  agreement. 


•).cn 


b)  Comparison  of  modal  and 
network  raaponsaa 


Fig.  4.  Second  Order  Modol  Identification  of 


a)  Input  -  Output 
eignala  for 


b)  Comparison  of  model  and 
network  responses 


l.A 


I.M 


Fig.  3.  First  Order  Model  Identification  of 
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2.3  Program  Description 


This  FORTRAN  IV  program  determines  a  linear  model  (transfer 
function)  of  a  network  from  recorded  laboratory  responses.  The  linear  model  is 

obtained  via  the  pencil-of- functions  method  discussed  in  Section  2.1.  The 
program  has  certain  features  which  are  discussed  below. 

Network  modeling  involves  the  determination  of  the  coefficients  a.  and  B. 

X  1 

of  a  rational  transfer  function  of  the  form 

6  +  6,s  +  .**  +  S  s 

H(s)  =  -2 - ^ - D_  (17) 

a  +a,  s  +  ...+as 
o  1  n 

such  that  the  output  of  this  model  to  a  given  input  will  approximate  the  actual 
network  output  to  the  same  input.  Equivalently^  in  discrete  time  ll7  ]  we  wish  to 
determine  the  coefficients  a.  and  b.  of  a  function  of  the  form 


H(2) 


b  +b,2^+...+b2 
o  1  n 


a  +a,2^+...+a2 
o  1  n 


(18) 


If  the  network  under  study  is  assumed  to  have  direct  transmission »  the 
numerator  coefficient  b^  is  nonrero.  This  choice  of  model  structure  is  imple¬ 
mented  by  setting  IREM  =  0.  I'Jhen  direct  transmission  cannot  be  assumed  (i.e., 
it  is  known  on  physical  grounds  that  the  impulse  response  of  the  network  will 
not  contain  an  impulse),  then  b^  should  be  set  to  2ero.  This  is  accomplished 
with  IREMj^O.  For  example,  if  IREM=1,  the  coefficient  b^  in  equation  (18)  is 
set  to  aero;  for  IREM=2  the  coefficients  b^  and  b^^  are  set  to  aero.  It  is 
recommended  to  use  1REM=1  whenever  direct  transmission  cannot  be  assumed. 

All  calculations  are  performed  in  discrete  time;  finally  H(z)  is  transformed  by 
means  of  a  pulse  invariant  transforraation*(IZTS=2)  to  the  corresponding 

continuous  time  model  H(s).  After  modeling  has  been  accomplished,  the 

li 


*  See  Appendix  D 


Normalized  mean-square  error  fend  its  square  root)  comparing  the  model  and  actual 
network  responses  are  calculated  (subroutine  ERROR).  These  errors  are  calculated 
as  shown  below. 

I  '■'«'>  - 

N.M.S.E.  =  - - - 

S  x^(k) 
k 


R. N.M.S.E.  =  .'N.M.S.E. 

Another  feature  of  the  program  is  the  capability  for  blas-remov.'.l  from 
the  recorded  laboratory  responses  (IB1AS=1).  This  feature  allows  consideration 
for  bias  that  may  have  been  introduced  through  the  laboratory  measurement 
system  to  the  recorded  output-input  data. 

Finally,  a  plot  option  (IPLT)  is  available,  l^hen  1PLT=1,  two  sets  of 
plots  are  given.  The  first  shows  the  original  output-input  data  measured 
from  the  network.  The  second  plot  contains  the  original  network  response  and 
the  identified  linear  model  response.  This  plot  allows  visual  inspection  of 
the  closeness  of  the  model  fit  to  the  actual  (desired)  response. 

To  enable  the  test  engineer  to  effectively  use  the  program,  a  description 
of  the  input  data  cards  is  given  below. 


I 


rS 


INPUT  DESCRIPTION 


CARD  #1  The  first  card  is  a  title  card.  Columns  1  through  80  are 

available  for  an  alpha-numeric  title. 

CARD  #2  Option  card  containing  three  variables 


Variable  Name 
(Format) 

Description 

Columns 

N(15) 

Order  of  the  system 

1-5 

MP1(15) 

Number  of  data  points 
(output-input  data) 

6-10 

IPLT(I5) 

Plotter  option; 

IPLT  =0  no  plots 

IPl.T  =  1  plots  on  line  printer 

11-15 

D\RD  through  CARD  12-fN0UT| 

NOUT  =  ((MPl+7)/S],  where  {X]  is  the  truncated 
value  of  X. 

The  output  data  is  entered  on  these  cards  in 
8F10.0  fields. 

CARD  g[3+N0UT]  through  CARD  |2+N0U I+NIN] 

NIN  ®  l{MPl+7)/8],  where  (X)  is  the  truncated 
value  of  X. 

The  input  data  is  entered  on  these  cards  in 
SFIO.O  fields. 


"CARD  #[3+N0UT+XINl  Second  option  card  containing  six  variables. 


Variable  Name 
(Format) 

Description 

Columns 

N(15) 

Order  of  the  system 

1-5 

MP1(I5) 

Number  of  data  points  (output- input  data) 

6-10 

ISK1P(15) 

This  variable  determines  the  sequence  of  points 
plotted  on  the  printer.  If  ISKIP  =  1  every 
data  point  is  plotted,  and  if  ISKIP  =  5  every 
fifth  point  is  plotted,  etc. 

11-15 

1RE.M(15) 

Variable  used  to  specify  model  structure  for  the 
identified  system.  If  direct  transmission  is 
assumed,  iREM=0.  For  IRDI=m,  the  first  m  terras 
(in  ascending  order)  of  the  model  numerator 
polynomial  are  set  to  zero.  It  is  recommended 

1REM=1  when  direct  transmission  cannot  be  assumed. 

16-20 

*At  first  glance,  this  card  may  seem  partially  repetitious  with  CARD  f2. 
However,  when  multiple  identification  runs  on  the  same  output-input  data 
are  desired,  then  more  than  one  such  option  card  may  be  placed  here,  with 
the  option  variables  changed  as  desired  (for  instance,  a  run  on  only  part 
of  the  output-input  sequence  may  at  tiroes  be  needed). 
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IBIAS(I5)  Bias-removal  option  21 

IBIAS  «  0  no  bias  is  assumed  present  on  the 
output-input  data. 

IBIAS  =  1  bias,  assumed  to  have  been  introduced 
by  the  measurement  system,  is 

removed  before  identification  is 
performed. 

DELTA(F5.0)  Sampling  interval  26 

END  OF  FILE  CARD 


A  listing  of  the  FORTRAN  programs  used  is  given  in  .Appendix  A  . 
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III.  COMPUTER  ROUTINE  FOR  HIGH  PRECISION  INVERSION  OF 
SECOND  ORDER  VOLTERRA  RESIDUE  MATRICES 


The  determination  of  the  residues  of  the  quadratic  TF,  H^iSj^.s-),  in¬ 
volves  the  solution  of  a  set  of  linear  equations.  Unfortunately,  the  nusber 
of  equations  involved  are  large,  for  example  12  [2]  even  for  a  modest  single 
pole-pair  situation  (l.e.,  where  the  linear  TF  has  two  distinct  poles). 
Solution  of  these  equations  can  lead  to  computational  errors  unless  extreme 
caution  is  exercised  in  the  Inversion  of  the  associated  matrix.  In  fact  the 
problem  is  further  aggrevated  in  cases  where  the  system  is  wide  band,  i.e., 
when  the  poles  of  H^(s)  are  spaced  several  decades  apart.  In  such  situations, 
the  poles  of  H2(s.,S2)  Involve  sums  and  differences  of  the  linear  TF  poles 
which  can  result  in  precariously  close  values.  For  example,  if 

=  50  radians/s 

=  50  M  radlans/s 

then 

Aj^  +  A^  =  50.00005  Mradlans/s 

Aj  -  A2  =  49.99995  Mradians/s 

This  in  turn  causes  the  associated  columns  of  the  coefficient  matrix  corre¬ 
sponding  to  these  poles  to  be  almost  scalar  multiples  of  each  other.  The 
matrix  thus  becomes  nearly  singular,  or  highly  ill-conditioned  to  invert. 

The  program  presented  in  this  section  is  designed  to  deal  with  such 
wideband  cases,  and  more  generally,  to  invert  ill-conditioned  matrices  where 
ever  they  may  arise.  It  is  hoped  that  by  mastering  the  various  capabilities 
of  this  routine  the  analyst  can  cope  with  almost  all  situations  of  practical 
interest. 

The  program  possesses  the  following  features  which  enable  high-precision 
inversion: 

Adaptive  Scaling 

Application  of  Perturbation  TTieory  to  Ill-Conditioned  Matrices 

Iterative  Correction 

Before  discussing  each  of  these  in  detail,  it  is  useful  to  define  the 
term  "ill-conditioned"  matrix  -  We  will  call  a  matrix  ill-conditioned 

if  (a)  the  rows  (or  columns)  of  the  matrix  are  nearly  dependent,  (b)  "small" 
changes  in  one  or  more  entries  of  the  matrix  result  in  large  changes  in  its 
inverse,  or  (c)  the  nonzero  entries  of  the  matrix  differ  widely  by  several 
orders  of  magnitude  (and  remain  so  even  after  appropriate  scaling  has  been 
performed) I  18) ,  119),  [20).  Note,  the  above  conditions  are  not  mutually 
exclusive. 
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3.1  ADAPTIVE  SCALING 


In  many  applications  the  entries  of  a  matrix  differ  widely  in  their  respec¬ 
tive  sizes.  For  a  linear  system  of  equations  this  situation  arises  when  the 
values  of  the  (unknown)  variables  are  orders  of  magnitude  different  and/or  the 
various  equations  have  right-hand-sides  which  are  orders  of  magnitude  different. 
This  situation  can  be  remedied  in  many  cases  as  follows.  Denote  the  matrix 

of  Interest  as  A  .  Then  it  is  possible  to  factorize  A  as 
o  o 

A  =  PA  =  PAQ  (1) 


o  ‘“1 

where  P  and  Q  are  suitable  diagonal  scaling  matrices  [19  ] .  The  following  method 
was  developed  to  obtain  the  diagonal  matrices  P  and  Q,  and  hence  the  new  matrix, 

A. 

The  diagonal  entries  of  matrix  P  are  successively  computed  from  the  product 
of  ail  "significant"  terms  In  the  successive  rows.  The  term  "significant"  can 
be  specified  by  the  user  (in  the  examples  presented  here  any  entry  greater  tluin 
15  orders  of  magnitude  below  the  largest  entry  in  the  row  of  interest  was  con¬ 
sidered  significant) ; a  default  value  of  15  orders  of  magnitude  is  assumed.  Then, 
the  P^j^th  entry  of  the  diagonal  matrix  P  is  computed  as  the  (n^)th  root  of  the 

magnitude  of  the  aforementioned  product,  where  n,  is  the  number  of  terms  in 
1  ^ 

the  product. 


I, 


The  scaling  of  tl\e  1th  row  may  be  stated  mathematically  as  follows;  let 

a  =  MaxABS(A  ).,  [largest  entry  of  1th  row) 
i  j  o  ij 

*  10  [threshold  for  ith  row  (m  choof-eu  by  the  user)) 

(A  ),,  :  ABS(A  ).,>S^  [qualifying  entries  of  ith  row) 
o  i.i  oil  1 


n^  =  number  of  qualifying  entries  in  1th  row 

Then  (  1 

Pm  =  ” 

qualifying  j 

entries  ^  ' 
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At  tills  point,  wo  luivo  factorized  matrix  A  Into  two  matrices,!  .o.,  A  = 
'  o  (> 

1’  *  Aj,  wlioro  tlio  entries  of  A^  are  specified  as  follows: 


‘■'I'lj  ■  <‘'.>i.i'''ii 


I . n 


It  should  he  noted  that  in  certain  cases,  the  above  row  scaling  will  suf" 
fice,  and  further  sealing  may  not  be  necessary.  However,  In  general,  the 
above  process  may  be  repeated,  this  time  utilizing  eoluiiiii  scaling.  Spe¬ 
cifically,  the  column  scaling  Involves  factorizing  A^  such  that  A^  “A  *  Q, 
where  Q  is  diagonal.  'I'he  entries  of  A  are  obtained  as 


1 . u 


The  entries  of  the  scaling  matrix  Q  are  chosen  In  the  same  manner  as  those  of 
P,  e.xccpt  that  columns  rather  than  rows  (of  A^)  arc  examined. 

Utilizing  this  technique,  the  desired  inverse  is  seen  to  be 

A"^  =»  q"^  a"’  P~^.  (4) 


Kxiunple  1 


0. 1 0000000  K+03 
0.19999989  H+Ob 
0.10000000  b+lO 


0.20000000  i:-04 
0.40000000  E-01 
0.10000000  E+03 


0.29999999  11-01 
0.60000000  11+02 


Inspection  of  matrix  A  shows  that  its  entries  differ  widelv  In  relative 

o 

magnitudes;  in  fact  there  is  a  difference  of  fourteen  orders  of  magnitude. 
Therefore,  scaling  can  bo  used,  yielding  the  following:  A^  «  PAj  “  PAQ,  where 
0.391486768E-01  0.0  0.0 

P  =  0.0  0.7829 7 339K+02  0.0 


Row  Scaled 
Matrix,  A, 


0.2554 364 7K+04 
0.2554363911+04 
-0.3162277611+04 


0.27427647E+04 


0.51087295K-03 

0.51087304E-03 

0.3162277611-03 


0.43538684E-03 


0. 3162277611+0 6j 

0.7663094311+00 

0.7663095611+00 


0.7663094911+00 


Row  and  column  scaled  matrix, 


I  '  ^4’ 


H  j.S -S I  f \ '  e-S®*! 


0.93131018E+00  0. 11733771E+01 

A  =  0.93130987E+00  0. ] 1733773E+01 

-0.11529525E+01  0. 726314477E+00 

Now,  =  Q  ^  A  ^  5*  and  the  inversions  yield: 

0 . 20000000E+05  -0 . lOOOOOOOE+02 

A  =  0.20000000E+12  -0. lOOOOOOOE+09 

o 

-0.19999996E+09  0. 10000000E-K)6 

The  product  A^(A^)  ^  is  given  by 

O.lOOOOOOOE+01  -0.17053025E-12 

A^(A^)“^  =  -0.95367431E-06  0. lOOOOOOOE+01 

_0.0  -0.95367431E-06 

Checking  Product  of  Computed  Inverse 

After  having  computed  the  inverse.  It  is  desirable  to  check  the  accuracy 
of  inversion.  The  obvious  way  to  do  so  is  to  obtain  tlie  product  of  A^  and 
and  to  examine  how  close  it  is  to  the  unit  matrix.  However,  in  view  of 
the  finite  word-length  of  the  computer,  this  product  must  be  computed  care¬ 
fully.  Whenever  scaling  methods  are  employed,  the  product  of  A^  and  its  in¬ 
verse  may  be  determined  in  a  number  of  ways.  This  will  be  discussed  for  tliree 

different  cases:  (i)  row  scaling,  (li)  column  scaling,  (ill)  both  row  and 
column  scaling. 

(i)  Row  Scaling;  A^  =  PA 

In  this  case,  A  ^  =  A  ^  P  Therefore,  as  a  check  on  the  dependability 
-1  °  -I 

of  A^  ,  the  product  A^A^  would  probably  be  examined  as  follows: 

A  A“^  =  P(A  A“^)  P~^  (5a) 

0  0 

on  the  other  hand, 

a"^A  =  a“^(P"^P)A  =  a"^A  (5b) 

O  O 

Equations  (5a)  and  (5b),  while  representing  the  same  quantity  (A^A^^  =  A^^A^) , 
may  not  be  found  equal  due  to  the  available  computer  accuracy  (finite  word- 
length).  Example  2  illustrates  this  point  with  an  extreme  case. 


().999999910E-H}0 

O.IOOOOOOIE+Ol 

0.0 


-0.54977761E-18  | 
0.99999999E-02  ' 
-O.66666666E-O5J 

0.26469779E-22  \ 
0.54210108E-19  | 
O.IOOOOOOOE+Ol  ! 
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EXAMPLE  2: 


A  = 
o 


O.lOOOOOOOE+03 

0.19999990E+19 

-0.50000000E-08 


0.20000000E-0A  0. 29999999E-01 
0.40000000E+12  0.60000000E+15 
0.49999999E-11  0.0 


Factorizing  A  ,  we  obtain  A 
°  o  o 


PA,  where 


P  = 


0.39148676E-01 

0.0 

0.0 


0.0 

0.78297339E+15 

0.0 


A  = 


0.2554 364 7E+04 
0.25543639E+04 
-0.31622776E+02 


0.51087295E-03 

0.51087304E-03 

0.31622776E-01 


Now,  matrix  A  is  inverted,  yielding 


0.78297352F.+03 

0.78297352E+06 

-0.26104324F.+07 


-0.78297339E+03 

-0.78297339E+06 

0.26104333E+07 


.-1 


The  product  AA  is  conipvited  as 


0.99999999E+00 

-0.23283064E-09 

-0.45475735E-11 


0.2 32 83064 E-09 
0. lOOOOOOOE+01 
0.27284841JE-n 


The  required  inverse,  A^\  is  computed  as  A^^ 


0.0 

0.0 

0.15811383E-09 


0.75630943E+00 

0.76630956E+00 

0.0 


-0.22097009E-14 

0.316227766E+02 

-0.210818510E-0] 


-0.43368086E-17 

-0.43368086K-17 

0.99999999E+00 


A"^P'^ 


-1 


Once  A  Inis  been 
o 


computed,  we  can  obtain  the  product  of  A  and  A  by  use  of  equation  (5a)  or 


o 


(5b).  Utilizing  Eq.  (5a),  yields 


m=7 


I 


Comparison  of  this  matrix  with  the  unit  matrix  (the  (2,3)  entry  in  partic¬ 
ular)  would  yield  the  faulty  conclusion  that  an  accurate  inverse  lias  not 
been  obtained.  If,  on  the  other  hand,  Eq  (5b)  is  used  as  a  check,  we  find 
that 

0.99999999E+00  -0.6987(>K5(iE-Ui  -0. 170530251;- 12 


:v"^A  =A“^(P~’'t>)A=A  ^A  = 
o  o 


-0. 2979532.'.  E-06 
0.2542S985E-05 


0.99999999E+00  -0. 101 8()340E-09 

0. 36878730E-12  0. lOOOOOOOE+01 


Comparison  of  this  matrix  with  the  unit  matrix  would  bo  favorable.  Tims,  in 

the  case  of  row  sea  ling,  equation  (5b)  must  be  used  .is  a  check  on  the  product 

A^^  and  A^^  due  to  the  finite  computer  word-length  .  Use  of  this  equation 

avoids  the  possible  problem  (present  in  I'.q.  5a  )  that  the  matrices  1',  (AA  *) 

and  P  ^  may  not  be  eompatible  for  multiplication,  as  exemplified  above.  To 

reiterate,  in  the  c;isc  of  row-scaling  the  product  of  A  and  A  '  must  he 
-1-1 

computed  as  A  A  =  A  A,  where  A  is  the  scale.l  matrix. 

^ ''•g:  lit  a  similar  manner  we  will  show  that  the  product 

AA  ^  must  be  used  as  a  cheek  on  the  goodness  of  A^^'  in  the  case  of  column 

scaling.  This  can  bo  seen  from  .i  eoinpar isoii  ,'f  tbe  lolli>wing  equations, 

where  A  =  AQ  and  A  ^  =  Q  *A. 
o  o 


a"^\  =  Q"‘(A"^\)t) 

o  o  ^ 

A  a"’  =  A(QQ”^)a“’ 

o  o 


AA 


(tia) 

((ibl 


Again  due  to  finite  word  lengtlt  in  the  eomputor,  equation  (bb)  should  bo 
used  to  verify  the  goodness  of  t  lie  inverse  when  eolutmi-seal.ug  is  performed. 


I'M 


si* 


i  - 

- 


ili)  Row  .Hid  Column  8_<.vi_l  i cig:-  in  this  case,  matrix  A  is  factorixed  as 
.  '  _] 

A  =  PAQ,  with  the  required  inverse,  A  =  Q  A  P  .  The  dee  is  ion  to  use 
-1  -1 

A  A  or  A  A  as  a  clieck  towan.  llto  .accur.iev  of  tlic  inverse  obtained  mav  he 
o  o  o  o  ■  _J 

arrived  at  as  follows.  The  equations  represent ing  A  .A  and  A  A  are 

2  3 


=  PA(QQ“^)A”^p“^  ■=  P(AA“^)P~^ 


a"^a^  =  Q~^a“^(P“^P)AQ  =  Q“^(a“\)Q 


The  difference  In  magnitude  between  the  largest  and  smallest  entries  in  tlie 
diagonal  matrices  P  and  Q  is  calculated.  These  two  differences  are  com¬ 
pared,  and  the  matrix  with  tlie  largest  difference  is  attempted  to  be 
eliminated  (either  P  or  Q) .  Thus,  if  the  entries  of  P  are  more  incomr,atible 
for  multiplication  tlian  tliose  cf  Q,  equation  (7b)  would  be  used  as  a  check 
rather  than  equation  (7a).  An  example  will  help  to  clarify  this  procedure. 

EXAMPLE  3: 


O.lOOOOOOOE+03 

0.20000000E-04 

0.29999999E-0r 

A  = 
o 

0.19999990E+19 

0.40000000E+12 

0.60000000E+15 

-0.50000000E-08 

0.49999999E-11 

0.0  i 

matrix  is 

now  expressed  r-s 

A^  =  PAQ,  where 

0.39148676E-01 

0 

0.0 

P  = 

0.0 

0.78297339E+15 

0.0 

0.0 

0.0 

0.15811388E-09 

1 

1  0.59091075E+03 

0.0 

•1 

0.0 

Q  = 

1  0.0 

0.20208866E-02 

0.0  ! 

1  0.0 

0.0 

0.76630949E+00_ 

1 

1 

0.43227589E+0] 

0.25279643E+00 

O.9999999IE+O0' 

A  =  1 

0.43227575E+01 

0.25279647E+00 

O.lOOOOOOOE+01 ' 

1 

-0.535153177E-01 

0.15647971E+02 

0.0  ! 
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-0.20003991E+07 


0.20003998E+07 


-0.16155222E-01 


Comparing  the  magnitude  difference  between  the  largest  and  srallest  entries 

■’4 

of  P  and  Q,  it  is  found  tl>at  in  P,  the  difference  is  approximately  10  , 

while  in  Q  the  difference  is  approximately  10^.  Therefore,  Eq.  (7h)  sliould 
be  used,  in  order  that  multiplication  involving  P  and  P  ^  is  avoided. 
Performing  tlie  required  multiplication  indicated  in  Eq.  (7h),  we  obtain 


a"^a  =q'^(a"‘a)q 

o  o 


0.99999999E+00 

-0.23651533E-06 

0.89538983E-06 


-0.70030334K-16 

0.99999999E400 

0.15661482E-12 


-0. 5661 39 84 K- 13  j 
-0.43109444E-10  j 
0.  lOOOOOOOE+01 ! 


3.2  PERTURBATION 

In  some  applications,  th-'  given  matrix  may  bo  iU-c«'udi tinned  to  the 
point  that  the  scaling  method  described  above  will  not  allow  inversion  of 
the  matrix  to  the  desired  precision.  In  this  case,  application  of  per¬ 
turbation  theory  to  th<'  sealed  matrix.  A,  may  be  helpful.  Several  methods 
will  be  described. 


A.  Diagonal  P^turbation 

The  first  method  consists  of  forming  a  new  matrix  (IS) 

C  *  A  +  CO  (81 

where  A  is  the  scaled  version  of  t'  >iiginal  matrix  and  1)  is  a  diagonal  matrix 


- T . . . . . . . . . . . . . . . . . . . . . . . . . . . 'yTOafiRr’' . 


whose  entries  can  be  taken  as  those  of  the  diagonal  of  A.  The  nuiltiplier,  ! , 
is  chosen  to  he  sviitahly  small  (more  will  ho  said  about  this  choice  later) 
such  that  C  is  invertible  (allowing  for  the  available  computer  accuracy).  We 
can  now  write  A  =  C  -  rD  so  that  we  have  parameterized  A  in  terms  of  the 
small  parameter,  r,.  Since  A  =  A(f)  is  an  analytic  function  of  r,  its  inverse 
is  also  an  analytic  function  of  c;  therefore,  a  power-series  may  be  written. 


Indeed , 


A~^  =  (0  -  CD)~^ 

=  (1  -  CC"^1))"V^ 


-1  -1  -1  ?  -1  ?  -1 
=  (C  +  CC  DC  +  G  (C  D)  C  +  ...] 


Thus,  using  C  ,  D,  and  c,  the  inverse  of  A  may  be  computed.  An  advantage 
of  this  method  is  that  it  can  be  performed  without  any  visual  inspection  of 
the  given  matri.\.  This  is  so  because  the  diagonal  matrix  D  is  specified 


automatically.  A  disadvantage  worth  noting 


is  that  the  required 


inverse,  A  is  represented  as  an  infinite  series  of  matrices  and  cannot  be 
expressed  in  closed  form.  Thus,  an  exact  representation  of  A  ^  cannot  be 
obtained,  althougti  it  can  be  approximated  to  the  required  accuracy  by  computing 
and  summing  enough  terms  of  the  series. 

We  now  return  to  the  matter  of  choosing  a  suitable  value  of  C.  Clearly, 
the  smaller  the  values  of  c,  the  faster  *>ill  be  the  convergence  of  the  series 
(9).  On  the  other  hand,  c  must  be  large  enough  so  that  it  results  in  adequate 
perturbation,  i.e.,  such  tliat  C  is  Invertable  with  the  available  computer 
accuracy.  It  is  usefvil  to  note  that  a  theoretical  upper  bound  on  c  can 
be  obtained;  Indeed  for  the  series  ( 9  )  to  converge  it  is  necessary  that  Ifj  be 
less  than  1.0/ [largest  eigenvalue  of  C  ^d|  .  119),  120). 

B.  .Single-entry  Perturbation 

An  alternate  perturbation  method  consists  again  of  forming  a  new  matrix 

C  =  A  +  CD 

where  we  now  restrict  the  perturbation  matrix  D  to  have 


only  one  non-soro  entry,  again  on  the  diagonal.  Without  loss  of  general iiy, 

this  non-zoro  entry  can  be  taken  to  bo  unity.  Then,  the  r.iatrix  1)  lias  rank-one, 

trace-one  and  hence  can  be  written  as 

D  >  XY^  (10> 

-1  T  -1 

A  »  (C  -  rXY*)  ^ 

-1  T  _1  -1 

-  |1  -  .  (C  ‘x)y‘  1  ‘c 

-  1  T  - 1  T  - 1  T  - 1 

1  +  rC  XU  +  rY  C  X  +  (rY  C  X)‘  +  . . .  hY  C  )  (IH 

Now,  the  quantity  in  brackets  in  eqxiation  (lO  is  a  power-series  involving 
scalar  terms.  Therefore,  equation  (11)  c.au  be  expressevi  as 


-1  r  1  "1  T  -1 

+  rC  ‘X  - ‘--y-  (y'c  ') 

_1  -  fY‘c~‘x_ 


It  can  ho  shown  that  the  series  converges  for  all  r.  It  is,  however, 

best  to  choose  f  — f  •  Thus,  equation  (12)  represents  an  e.xact  solution 

Y  C"^X 

for  A  .  A  disadvantage  of  tltis  irothod  is  that  visual  inspection  of  the  matrix, 
A,  is  necessary  to  determine  the  row(s)  and/or  column(s)  causing  difficultv  in 
the  inversion  process,  and  thus  a  suitable  entry  to  perturb.  An  illustrative 
example  of  this  method  is  given  below. 

KXAMl’LK  4: 


O.lOQOOOOOK+03 
0.19999989E+06 
-0.  IOOUOOOOF,+10 


0. 20000000E-OA 
0.40000000K-01 
0.10000000E-K)3 


0.29999999K-01 

0.60000000E+02 


A  visual  Inspection  of  this  matrix  shows  that  tltc  first  and  second  rows 
are  nearly  identical  to  within  a  multiplicative  factor  of  ".0  x  lO"^.  Therefore, 
suit.able  diagonal  entries  for  perturbation  would  be  either  the  (1,1)  or  (2,2) 
entries.  Tlie  (2,2)  entry  was  selected  for  perturbation,  and  through  iteration, 
a  suitable  value  of  C  was  found  to  be  r  »  1.  x  10  Scaling  was  first  per¬ 

formed,  and  then  the  portnrbatioiv  yielding; 


27 


0.931.3I018E+00  0. 11733771E+01  0.99999991E+00|  0.0  0.0  0.0 

A=C-eD  =  0.93130987E-K)0  0. 11733773E+01  0. lOOOOOOOE+01  -10^°  0.0  1.173377  0.0 

-0.11529525E+01  0. 72631447E+00  0.0  0.0  0.0  0.0 


T 

Note  that  0  can  be  expressed  as  D  =  XY  so  that 


D  =  11.173377 


1  [o. 


0  1.0 


Tlie  inverse  of  the  original  matrix  can  now  be  computed  as  discussed: 


0.20000000E+05 

A~^  =  0.20000000E+12 
o 

-0.19999996E+09 


-O.lOOOOOOOE+02 

-O.lOOOOOOOE+09 

0.10000000E+06 


0.52504835E-18 

O.IOOOOOOOE-Ol 

-0.66666666E-0S 


This  yields  the  product 


"o.i^ 

A  A~^  =  -0.9 

o  o 

0.0 


O.lOOOOOOOE+01 


0.95367431E-06 


-0.51159076E-12 
0. 10000000 E+01 
-0.95367431E-06 


0.66174449E-22 


0.99999999E+00 


C-  perturbation:  THE  LIMITING  CASE:  In  the  previous  discussions  on 
perturbation,  it  was  shown  that  matrix  A  is  a  function  of  matrix  C  and  the 
scalar  quantity,  r..  That  is  A  =  (C(e),c).  Recall  also  that  wliereas  A 
was  ill-conditioned  (for  inversion),  there  were  certain  values  of  c  for  which 
the  newly  formed  matrix  C  could  be  made  well-behaved.  Inspection  of  equation 
( 8 )  shows  that 

A^=£imC^  ft -w 


t  I 

I  I 


This  observation  can  be  exploited  in  the  following  way.  Successively 


small  values  of  €,  say  e^,  are  used 


to  form  a  family  of  C  “  C(c^) 


matrices.  The  inverse  of  C(e^)  is  computed  for  each  value  of  used,  and 
the  successive  Inverses  are  examined.  There  will  exist  a  region  wherein 
reducing  the  value  of  e  from  to  will  have  little  effect  on  the  entries 
of  C  This  is  shown  grapliically  in  Figure  1. 


.-1 


(Shaded  region  is  where  C 
is  not  well-behaved) 


^min  ^i 


small  r 


large  C 

KIg.  1.  Effect  of  small  changes  in  e  on  C 
In  this  region,  C  ^  can  be  taken  as  an  approximation  to  A 


Figure  1,  there  will  exist  some  t  . 

®  min 

behaved.  The  closer  the  selected  value  of  •'  is  to  this  value  of  f 

‘  ^  Although  this  method  may  only 

yield  an  approximation  to  the  actual  required  inverse.  A  \ 


As  seen  in 

for  wl>ich  the  inverse  of  C  is  well- 

Ihe 

_ j  min 

better  will  become  the  approximation  A  =  C 

it  may  bo  further 

refined  by  use  of  the  method  described  in  the  next  section.  In  fact,  the 
iterative  correction  method  (of  the  next  section)  may  be  used  in  conjuction 
with  all  of  the  methods  discussed  earlier. 

3 . 3  ITEfcXTlVK  CORRECTION 

Consider  a  matrix,  X,  and  assume  that  its  inverse  has  been  computed  as 
Y  X~^  .  The  iterative  eorrectior  method  (see  Fig.  2)  consists  of  forming  the  product 
XY,  comparing  it  with  the  identity  matrix,  and  improving  the  computed  inverse, 

Y,  by  an  amount  proportional  to  the  error  between  XY  and  the  unit  matrix.  To 
examine  the  effect  of  this  operation,  let 

(lA) 


Y  »  X“^(l  +  E) 


where  1  =  identity  matrix  and  K  is  equal  to  the  difference  matrix  between 
XY  and  1 . 

2*} 


XY  =  (I  +  E) 


(15) 


E  =  (XY  -  1) 

Now  consider  the  iteration  Y.  =Y  -  YE  (121.  Cle^lrlv 

improved 

Y.  .  =  Y  -  YE 

improved 

=  X"^(l  +  E)(l  -  E) 

-1  ■> 

=  X  -  K-) 


(16) 

(17n) 


(17b) 


Upon  the  second  iteration, 
Y 


=  X~^(l  +  E^) 


improved 

and  so  on.  This  prodecure  can  be  depicted  in  block-diagram  form  as  in 
Figure  2. 


-E 


=c> 


MULTIPLIER 


r— 


XY 


MULTIPLIER  ^ 


=>'• 


Fig.  2;  Block  diagram  representation  of  iterative  correction  method. 

The  number  of  iterations  rc  be  used  may  be  specified  by  the  user.  For 

this  work,  n  iterations  have  usually  been  used,  where  n  denotes  the  dimension 

of  the  matrix  in  question.  Note  that  a  more  general  version  of  (17a)  is 

Y.  ,  =  Y  -  S  *  YE  where  B  is  a  suitable  positive  fraction, 

improved 

EXAMPLE  5;  (Effect  of  iterative  correction) 


A  = 


0.10000000E403 

0.19999999E-H)6 

-O.lOOOOOOOE+10 


0.20000000E-OA 

0.40000000E-01 

0.10000000E-K)3 


0.29999999E-01 

0.60000000E^2 

0.0 


2k 

t*- 
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Assume  that  the  inverse  matrix  has  been  computed  as 


A-'  = 

0 

0.19999999E+06 

-0.99999998E+02 

-0. 1005869 3E-16  1 

0.19999999E+13 

-0.999999988409 

0.99999998E-02  | 

-0.19999999E+10 

0.99999998E406 

-0.66666665E-05J 

so  that 

A  A-' 

0.99999999E+00 

0.36379788E-11 

0.52939549E-22 

=  ■  -0.15258789E-04 

0.99999999E+00 

0.0 

o  o 

-0.11718750E-01 

0-15258789E-04 

0.99999999E+00  1 

j 

Now,  if  1  iteration  of  the  correction  method  is  performed. 


A 

o  o 


O.lOOOOOOOE+01 

0.0 

-0.1.'1625000E-01 


0.0 

O.lOOOOOOOE+01 

0.0 


0.0 

-0.54210108E-19 

0.99999999E-H)0 


with  N(=3)  iterations  performed. 


A  A  ^ 
o  o 


0.99999999E-K)0 

0.0 

-0.78125000E-02 


0.0 

0.10000000E401 

0.0 


-0.13234889E-22 

0.54210108E-19 

0.99999999E-K)0J 


It  can  be  seen  that  the  product  A^A^  is  approaching  the  unit  matrix.  The 
worst  entry,  (3,1),  in  the  product  has  been  reduced  to  about  60%  of  its 
original  value  in  the  3  iterations. 

The  improved  inverse  has  been  computed  as 


i-l 


0.19999999E406 

0.19999999E+13 

-0.19999999E+J.0 


-0.99999999E+02 

-0.99999999E+09 

0.99999999E+06 


-0.76657257E-17 

0.99999999E-02 

-0.66666665E-05 
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This  example  has  demonstrated  the  usefulness  of  the  correction  method  in  im¬ 
proving  the  calculated  inverse. 

3.4  APPLICATION 

Determination  of  the  kernels  of  the  volterra  model  can  be  accomplished  by 
identifying  the  various  link  transfer  functions  H^fs^),  K,(Sj^,s,),  etc.  The 
first  step,  of  course,  is  to  determine  H^fs^)  for  the  small  signal  linear  case. 
Next,  identification  of  H.j(Sj,s,)  is  attempted.  Hovever,  it  is  shown  in  [l] 
that  the  poles  of  H_(s,,s-)  are  given  as  y..  =  a.  +  X,  where  X.,  X,  are  the 

^  i  ^  IK  IK  Ik 

poles  of  the  linear  transfer  function  Hj(Sj).  The  residues  at  these  poles 
can  then  be  determined  by  solving  a  set  of  linear  equations  from  data  obtained 
for  larger  amplitudes  where  quadratic  effects  become  nonnegligible. 

The  set  of  equations  which  must  be  solved  (for  computing  residues  of 
the  response)  may  form  a  nearly  dependent  set  if  the  network  is  a  wideband  (i.e., 
its  poles  are  separated  by  two  or  more  orders  of  magnitude)  network.  Accurate 
inversion  of  the  corresponding  matrix  can  be  quite  a  formidable  task  in  this 


W-  H 


In  this  section,  a  12  x  12  matrix  generated  in  the  analysis  of  a  wideband 
network  is  examined.  The  equation  of  interest  is 

R  =  (A  )'H’ 
o 

where  R  =  column  vector  of  residues  of  poles  of  the  system  response 

Y  =  column  vector  of  integrated  system  outputs  at  time  T 

A  =  12  X  12  matrix  whose  entries  are  generated  bv 
o 


(A  )..  = 
o  ij 


i  m-1 

I  {  - - ) - - - 

m=I  (m-1)!  X. 

1 


Tlie  linear  portion  of  the  system  examined  has  two  poles.  These  poles  are 
Xj  =  -0.01l550998(2r)(l0^)  rad/sec 

X,  =  -10.616986  (2-) (10^)  rad/sec 

The  system  is  excited  by  the  input  function 

X.  (t)  =  le“l^  +  u(t) 

in 


where 


-10  rad/sec 


-1.75  X  10  rad/sec 
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3.i7;o.327;j;i45i2:7n*33 
0.7G5i./.3w)2j2j;G.)P-32 
o.iCG.'vGG  )3i::i>3^::vD*3i 


0.7495JD410G?3J5J5S:d-j: 
0.35JD534  3‘»552il0Jjj3-31 
D.>G3C3jJ32Gl&3u772?0“Jl 


0.^4..1002GC02  273J5nu-32 
0.>030^G35707k<72?3o/k5*01 
0420'«35G37C3G<>303^^30*01 


0.35C31472G7152i2:4>r-31 

U*2S04lJoS7'j4gvG343/l>*37 

0.20544AU2I37Gyl327:u-02 


3.2DCoGW1;202C5{*>3:G-  *2 
0.12':;.'G7.»J2!.;i2?Cl%70-.* 
3,7..2;v0j231J734:55OC'-3^ 


0.>522J...>352*»5.1,)73Sr-)l 
0  »22HG/04Ji  /  rf2.»w274>D“»  I 
0*4«4/g4v4C30«3^4033gD*32 


0. 2  3i>»6G0720.f^&!>‘J^3D*02 
0.,>GJJi3;2330/ Jl;7^JD-32 


0.&3533J93:03520:C11/'02 

04&0yo/o50,GwJ3r27:0n-)3 

0.4051114713035321250-33 


0.50n33?05C:5155:Clf>.G3 
4.;2o:ui/53;53:c:i2d-v? 
J.2  a:7500524013271GD-32 


D.53072:G7:3«.J410;3:r-02 

0.44':0::47722350C:5:D-03 

0.:.0J.3g043315355;:7U-3> 


0.25  3,3  3211  u  722331073  D*  03 
D.15G5:21511430?5‘wv.')O2 
0.U3Gl?i354lj;52132H-02 


O.G4332Ci. U1223J4:7n-03 
0.735l303;j;3G5C;5‘.4I'-:4 
0. 5033352  537230330340-04 


.>.734  5024131702334:0')-:; 
u.  303740024 131 0007:50-33 
3, 1Q30G24 '3:370333:430-05 


0. 33:733700221  34G:G3D-.)5 
0.341.53171 7443:54-.>5'3-34 
3. 120/234101 3301 :334,'>- 35 


0.3052C1020051 34505^0-04 
0.2131310105/J5203.1D-05 
J.:5f4:d35075030.12:20-J5 


0.G43  :3009g00302;45G!5-34 
o.sf.047r.i4v0>div:o5if>3: 
0.C09051G:14027:7233D-35 


0.0G133C0100S5703012O-35 
0.:5o:J752u222C13435n-34 
0. 2244007  503G'i03319in- 34 


D.:40C.>4051C:C.l  007310- 34 
0,7G',4:374)‘.73093vh50-.*^ 
0.14:3070253.4103:400-34 


0.45C0132203732C3u. 00-35 
244530/5105433735:0-04 
G.17:752/45173413h-:0-34 


0,5520545335024397710-05 
0.041  4C4D91521552022I)-0G 
0.0065040031 2107G02u!'-00 


0.C4227C4013.712Ci30O-3G 
0. 3  J.;3?32  04220  7030460- 33 
0.211705G245193C3413D-05 


0,3 4 5>4 34 743050472300-33 
3.743G:j2G:c^»»^07373D'Ww 
0.1427;277..:;511 17000-;^ 


0. 4514557 lG25  3G3u7 330- DC 
0.2295374355960774.13-35 
D.170G3536.f54J3:75.5D-05 


0.33232144C737403437I5-3G 
0,7057JJi.3-:G23./55:73-:7 
0.57v0;  :5CilG4S33'3’.JI-37 


0./0035C05534?S07:.7rO-07 
0.  ^.i47v2533GyC573413I*-vv 
3.171:31473454294035j»-OC 


9.4125545230374013073-36 

O.G2G.7>/l..J3C*3i53lD-37 

0.1179375i2r.:7O27<.23-JG 


3.3:247202G3730l37jGn-07 

:.1v54>72732j:4j<^7:/D'0j 

0.133'324v02IjgJ37w^.^J-0o 


0. 4703271013034  37')53R-3: 
0. 51^11203  3:5226')454|- 03 
0.42C34-i52:522:5250f)-03 


O,51C570413232533G’5r- 
o.i7::3005io27:.co:/ir-o7 
0.122500  r755.'J5C72:2P-37 


o.27:i.'‘33iio2v>6;46in-:7 
9.40153G5Di.:^3iG13  )..{)-JG 
0 .G55325473'*G433 . 7 3.1'D 


0.2053C3*<o74*.301t3>3<.'-0v 
9.1515311u603.«3>  i«j4.>-j/ 
:.10J39i773227332jJ4C-O7 


0.G0795?5:52u53::C'.0n-05 
0. 3302744256^5 /sOG04n-03 
0.2791'39J1331G::1020!j-0J 


0.35:572>4510G7G3546D-09 
0.10j3J29:2  3)4  34.,4200-3G 
3.77:7.’15415553156:4D-39 


-9.20^^373)  4 3>37C3^C0U- ;3 
O.'.Jl.,  25r5G3/G;2:3U-39 
n.543o.2103953C56373U-03 


O.UG39129C4:3G<.11.3D-03 
^.03U:2J737  j/Jj34j3U-o  3 
0.  G4435G40  30432 ’.6223U-J3 


O.C3:C3C070C9:nG53D'»-94 

0.1‘3G3G0614243'J15:34'’i-10 

0.104720531/220:97551-10 


0. 1 030322  79?2:5::n9D-39 
3.5i::i‘33  35G24/4?579;»-13 
0.444  7G2  D0G75  35010C:»i}-l  J 


0.13149  33:5  K321375GU-07 
o.i;7:Gw:i1v^45»::/21)-io 
0,319  /353)6i:»JN:j}72n-lD 


0. 110:34367011133. )20-lj 
3.474:32:23044Oj11>3D-13 
•'«3/47l9/473./0<i57..jG()“2v 


-0.1154)  371003'J9303G4:)-)2 
o.ioc32ioi7ir.7/:jc)5:j-ii 
0.  <4  i;3  / 1502300122  7  Jo  3n- 12 


D.10:4:093/)Gg5222'»H)-i: 
0,307:29G44G0/0u3141O-n 
0.2  32  /505l5jG4777913n-i) 


-D.C9347G343:3:: 2:3940-97 
0.3579:j;5i:o53v.l2JD-l2 
i).16%279..327)  )2^::5n-ii 


9.53C124530223i:2:ilO-l2 
0.2j/;j434u7j223'».::50-12 
0. 1 3533213 4 5o7v4:<7  50-3 1 


It  is  readily  apparent,  through  visual  inspection,  that  row  scaling  is 
necessary,  yielding 


f  f  '  ( 
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(no;:)  scALCii  (UTiii;;  a 


0.1AlG431)«0ilA?:72tl*0: 

0.30125A12'335i<03'JG.'bliO0 

O.^GCUgOAJuJlA'JGOuOZIl'OO 


0. 03092 )00&03G51JJ:30*01 
a.GlA2G2G09ib3G)091GI>*d0 
0.330252»9a2b057A2S2li*00 


a.G2u'i7212349G3190GAIW31 
0.Gf-lOv9J2i79772Gb25n*JO 
0. 3C14i4 9012170 J97/onta3 


0.113&1173G5G01G370Gr)«01 

a,G003uS0G423u7G370GII*3O 

0.30G7SGG501bGg7792GiWU3 


O.GG20371940i)51390S31HOl 
a.G3S0OS22Sl'(335G037R'J0 
0.G0GGaG9330GG104  3SbI)*09 


0.393973397132595S920«01 

0.52G4217201G50SG032'i*j3 

0.4270uG54G30C177054l)*30 


0.35013SOG2'i337355:2n*31 
O.SGS)i317G23  3J2207(i'jIl«33 
O.GGG33ZU9050G500117r>  30 


0.30031220  739332004411*31 
0.SGu33071G9445794  2.ili*33 
3.4C30947u934CG392G7n*33 


0.4140203395034053420*32 

0.45G09115G37G459;C4!)*3a 

0.3753442331391970370*00 


-0.0200233374307023320*04 
0.3433437373301523310*03 
0. 2375103130513397530*  .)0 


3.5359511529057197190*00 

0.1272332137939*43340*00 

0.1353341*00332722020*00 


0.5359033152242729440*03 

0.539743317343)0199711-01 

0.4480722331333154320-01 


0.3010:.710*7JC.72iO:.33*09 
0.233  )424/2:.:i73)9)4n<O) 
0.123023:391202413270*01 


0.414CSj/177:.29149i3O*30 
0.2351754935749047040*01 
0.150030o3;.050:01  3490*31 


0.4515502121109233230*00 

0.227509:017449223:30*3: 

0.1541379557042:47210*91 


0.4:034  95..  71 351139220*00 
3.21971)-:i330944043D*01 
3.13:ii:9304403570:7n«01 


0.15::2527  )55124'.d;  ,D*02 
0.5145  )14, -0149211.  9.9*00 
0. *.92:5:1*2-75039*70*30 


9.::*9*7:?i.44  91515>;d*oi 
o.5:i32..9i.-5.;3:4:33in*oo 
0.7/201*.’14*u2?1.9;JO<OJ 


0. *197547145591 '.*4. ID*  01 
0. 39525590 J5CJJ:9J. -0*90 
0.:22Cu5] 709* 10:471 ID* 00 


0.53919*i3b7*71.!529JD*0) 
3.4221072*27*JJ****:D*93 
o.:s:g59* 9:4*71:11 . id*  do 


9. 1:9:255:73555* '4*123**0 
0.1377:5*?9120*29JJ4;i*j1 
0. *772* 9 9147**9* 97**0*99 


9.20 9997*07:1755995****9 
9.14*i412**3477**2J7?*91 
9.  *09 109* IS 03599* 3**. '*9* 


9. 7314:5*97  9379..  19  993*  00 
0.149257:1971*0121110**1 
).19205544410770l;0:0**l 


0.;49470;3555549/2o:d*90 

1).15u7295:7***935*4SD*91 

0.1053071750*9*1259*3*91 


0.5055925474340471050*0: 

0.21752:45:2552790;i:)*01 

0.1370253020957573410*01 


0.5:09732320035105370*00 

0.2000*14065159932430*31 

0.1573253771093344)40*31 


9.5459140333943917170*3* 

0.2a''243594S59:50339D*01 

3.137254139045153)900*31 


3.50335105350:3390730*33 

0.19*3732170570213900*91 

0.13/094430922)027350*91 


0. 450501 7374::i01374n* 05 
3.15323*2169222135500*01 
a.l07721C12C21C34932S*31 


0.34:05979215430433*0*33 

0.1094750533953771760*91 

0.7990421132134452:70*0: 


0.12734:9547597403590*39 
0.3:31401513525C1*93D*  9) 
0.2C45*757*025::35:an*0» 


9.540:3940023557*5240-9! 

0.150:713135015094793*09 

0.117;:i373103052.>540*01 


0. 45:539*150397 03:c:D* 91 
O.4'.5O9!9*J3*9.'4945-D*90 
0.0;.'.013-;751197*3.in*09 


0.30155399734949637*0*01 
9. 4:523513' !.’ 953:7 9. 0*03 
0.907C77934***C95552I)*30 


•1.5551354*31791072,70*31 
0.4:*2*77:0- 31 15*3:20*30 
0.92525745*47512C5):D*99 


0.55105:05*95)9031.10*01 

0.50590-597*992*97:20*99 

0.945295****911013**0*00 


0.2440277-:5931770. 50*31 
0.13011179*3 1 7 9;74-*D* 00 
0.7511 770C*1-.521213D* 00 


-0.2715951763:29442*10*30 
0.5!  i;..*7:J-079402910*09 
0.5;f.325-51*:/235145D*30 


0.341524349;*34174-3n*02 
0.1153'<>5:73  .077,2:90*00 
0.294) 595595 921 724 520*03 


-3.4535779).'?1493593)n*0'l 
0.4.'37.  5112:;:  9320490-01 
9. 55554  9594  9*70  9', *790- 31 


3.2o510*:14C35**194*D*99 

0.15019532931*!70*e.7.-'**l 

3.197079440:00290:5.'J*31 


3.279477317C:i7:47750*93 

3.149*35777)019**41*3*91 

*.19355702574949915*0*31 


0.29200739855**049:73*59 

9.143797>295S19412*-D*31 

0.1195959423492852270*91 


0.30:955509225*734410* 99 
3.1489159051119513**0*91 
0.1122:95927542:*4*j')*Jl 


3. 249021 13::94221 9*70* 39 
9.1157*0150029525977:1*91 
9.  .855054254951329:>D*93 


0. 1924040739*7 140::53* 39 
0.:5uC,4;93204958411r.*.9 
0.:G3227954*2*59475*'* 99 


0.7094733370:71215990-91 

3.3*3:152220359399210*93 

3.23733974510*57*9950*39 


9.505640726:01*4:5590-01 
0.125392359G2409:9*:D*90 
0. 991:130302:7*3471, '.-Ol 
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Program  Steps  Used: 

The  row  scaling  indicated  on  the  earlier  pages  was  achieved 
by  choosing  (see  program  description  on  page  39) 

ISP  »  1 

ISQ  =  0 

As  mentioned  earlier,  the  second  and  fifth  columns  of  the 
original  matrix  form  a  linearly  dependent  set.  Therefore,  a 
perturbation  procedure  is  needed  to  obtain  the  inverse.  We 
used 

XPERT  -  2 

so  that  the  2,2  entry  of  the  scaled  matrix  A  was  perturbed.  The 
amount  of  perturbation  used  was  FAC  =  O.lE-05.  The  inverse  of  the 
perturbed  matrix  C  was  obtained  by  the  high  precision  routine 
GKRDCT.  From  C  an  estimate  of  A  ^  was  obtained  by  use  of  the 
deperturbation  routine  DPERTl.  This  approximate  inverse  was  then 
refined  by  use  of  the  routine  C0RCT2  via  the  option  parameter 

ICORCT  =  1 

The  resulting  inverse  matrix  A  ^  was  then  descaled.  The 
inverse  of  the  original  matrix,  A^  ^  ,  as  well  as  the  product 
A^  ^A^  are  printed  on  the  succeeding  pages. 

(Note  that  a  preconditioning  of  the  original  matrix  was 
employed  by  setting  ICOND  =»  1.  This  caused  dia  /'.nal  entries  of 
A^  to  be  multiplied  by  1.000000001  »  1  +  l.OE-O^., 
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3. 5  Program  Description 

This  FORTRAN  IV  program  performs  high-precision  matrix 
inversion.  This  computer  program  has  cap.>bllitics  for  automatic  adaptive 
scaling  of  tlie  original  matrix.application  of  perturbation  techniques  in 
finding  the  inverse  and  iterative  correction  of  an  (approximate)  inverse  matrix; 
each  of  these  has  been  discussed  in  the  earlier  sections. 

To  enable  the  test  engineer  to  effectively  use  the  program,  a  description 
of  tlie  input  data  cards  is  given  below. 


Input  Data  Deck 

The  input  data  deck  consists  of  one  card  containing  input  variables, 
N^+T 

followed  by  (  — ^  1  cards  containing  the  entries  of  the  matrix  to  be  in¬ 
verted,  where  N  is  the  dimension  of  the  matrix  and  [,\1  is  the  truncation 
function. 


CARD  t'l  Option  card  which  contains  eight  variables. 

Variable  Name  Description  Columns 

(Format) 

N(I2)  dimension  of  the  square  matrix  to  be  inverted  1-2 

ISP(I2)  adaptive  scaling  option  3-4 

ISP  =  0  row  scaling  is  not  performed 
ISP  =  1  row  scaling  is  performed 

ISQ(12)  adaptive  scaling  option  5-6 

ISQ  =  0  column  scaling  is  not  performed 
ISQ  =  1  column  scaling  is  performed 

IPERT(I2)  diagonal  perturbation  option  7-S 

IPERT  =  0  perturbation  is  not  employed 


IPERT  =  1,2,. ..  ,N. perturbation  of  the  xpert 

entry  is  performed 

IPERT  =  N  +  1  perturbation  of  the  entire  diagonal 
is  carried  out 

■ ISLID(I2)  de-perturbation  option  to  be  used  with  IPERT  =  N+1  9-10 

ISLID=  0  deperturbation  is  performed  with 
subroutine  DPERT2 

ISLID=  1  deperturbation  is  performed  through  a 
"sliding'  correction  method  (utilizing  a 
family  of  matrices) 
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Columns 


Variable  Name  Description 
(Format) 

1C0ND(I2)  Pre-conditioning  option  11-12 

ICOND  =  0  r.:.  conditioning  of  the  original  matrix 
is  employed 

ICOND  =  1  the  original  matrix  is  preconditioned 
with  a  multiplier  1  +  1.0  E-9  along  the 
diagonal 

1C0RCT(I2)  option  to  be  used  in  conjunction  with  IPERT  =  0  13-14 

ICORCT  =  0  iterative  correction  is  not  used  when 
IPERT  =  0 

ICORCT  =  1  Iterative  correction  is  used  when 
IPERT  =  0 

lPRiNT(l2)  printing  option  15-16 

IPRINT  =  0  suppresses  printing  of  many  of  the 

intermediate  matrix  quantities  used  for 
computat ion 

IPRINT  =  1  all  intermediate  matrices  are  printed 

CARD  #2  through  card  ll(  ^  +  1] 

lliese  cards  contain  the  matrix  values,  and  are  read 
in  3D25. 18  format.  The  matrix  entries  are  entered 
on  the  cards  in  an  order  prescribed  by  columns.  That 
is,  they  are  entered  as  first  row,  columns  1  through 
N;  second  row,  columns  1  through  N;  etc. 

END  OF  FILE  CARD 


Note  a  logical  flow  diagr.am  depicting  the  effect  of  various  option  para¬ 
meters  is  given  below. 

A  listing  of  the  complete  FORTRAN  program  is  given  in  Appendix  B. 


PFRYIIRDATION 

(IP£RT>0) 

PERTURBATION 

(IPERI'O) 

•  SINGLE-ENTRY 

PERTURF'AT  ION 

DIAGONAL 

(0<IPER1<N) 

PERTURFATION 

(IPERT'-N)'s^ 

SINGLE-ENTRY  SLIDING 

DE-PERTURHATION  CORRECTION 
(OPERTI)  (IFERT'N, 

I  I  ISLID-I) 


01  AGONAL 

OC-PERTURRAUON 

(0PERT2) 


FINAL  BID 
FOR  CORRECTION 
(CORCl  I  OR  2) 


DE-SCALING 


NO  DL- 
PERTORBATION 
(ip;i.T-o. 
ICORCT=l) 


NO  DE-PFRTURBAl ION 
(IPIRT»0. 
IC0RCT”0) 


Flow  diagram  of  HPMINV  program  options 
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APPENDIX  A 


IGRAM  PROGRAM 
LISTING 


PROGfiA!;  ;JAKE;  IGRA.'l 

!!£TliORK  Tl:A!iS?ER-PU:iCTIO:i  in:«TIFICATIO:!  RCUTI'JE  UTILIZING 

pc:jcil-of-fu!ictio!is  tiETimn  (purc  imtcgrators  user) 

ni!‘.EIISIO:i  X(  1024 ),V(102A),XORG(lO2l!),yORG(  1024 ),xr.cc(102<i) 
Dlncnsio;!  DAlA(]02!i,2),RATA2(1024,2),r.UFF(  5372) 

OIJ!£!iSIO:i  r,(20,23),2(2O,70),CA(:J!A(2O),XI.Ai:nA(2!)),COrFF(2n) 
ni!;£i!Sio:i  liPutsscGii) 

0I(!£HSI0:J  TITLC(CO) 

REAL*C  G,2,G.MU-A,XtArtnA,C0:FF 
RC AL •  0  DC LTA, (1/  ".SAV,  PC LSAV, AVGQ, AVGW,  SUnV2 ,  XSA V 
EQUIVAtr.'lCE  (Z(l,l),RUFF(l)),(G(l,J),!!lirF(250a)5 
EOUIVALCriCE  (PATA(l.l),X(l)),(rATA(l,2),y(l)) 

EQUIVALCKCE  (XORGC 1  ),nATA2(l,l) ), <'.r;EC(  1  ),r)MA2 ( 1 , 2 > ) 
C0I"!0'1  /GKRD/IGKR 

coiHiou  /(Iu.".eu/:jn 

corj'.ou  /r.iAS/iGiAS 

RnCL=3.01 

|•.AX=20 

MAXPL=102li 

t)PT=HAXi’L 

URITE(G,2) 

'./l!ITn(G,102?) 

REAnC5,1021)TITl  E 

•.miTr(G,1021)TI7LE 

1;RI7E(C,1325) 

RE  AO  C ;  ,  1 00 1 ) .’-.P 1 ,  I  PL7 

RCAn(5,0095)(xor;r.(!;),K=i,.';Pi),(vnRG(i;),K=i,i:pi) 
p.r  AO  ( 0  ,  l ,  E-40^  1 2  5<I  > I!P1 , 1  SK I P,  I  RE.-S,  I  r.  I  AS,  nEl7  A 
IF(IJ.EQ.13700)GO  70  4320 
WR 1 77. C G ,  1 030  ):i,i;Pl , P7U7A,  I  REti,  I B I  AS 

igi;r=i 

IZ7S=1 

QSAV=1.0n00 

a=QSAV 

mii=:j-i 

I1P1=N+1 

NP2=!I*2 

np:ipi=:)*h*i 

Nr!iP2=';*t>2 

II»=tl-lREM 

ge:.’ERA7i:ig  sequeucc  x{k) 

t!02lil=l,MP] 
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2k 


6053 

C 

C 

C 


651*!. 

C 


100 

1231* 

C 

c 

1 

2 

66 

1003 


1001 

1003 

1004 
1021 
1022 
1023 

6905 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

p 

c 

c 

c 

c 

c 

c 


V(l)«VORG(t) 

x(t)=xop.G(n 

IrdFLT.EQ.OlGO  TO  0633 
WRlTn(6,1003) 

CALL  PLOTITCOATA  ,2,Mri,l,MPl,  I  SKI  P,fU4XPL,l,l  .0) 

COHTIIJUF. 

START  inE'ITIFICATIO:.’  FROM  Pipin’  OUTPUT  DATA 

CALL  GnAKI(X,V,K?l,'I.nELTA,Q^IZTS.GAia‘.\,XL.V.PA,G,Z,MAX,  IRF.’I) 
CALL  ERROR(XREC,yORf.,GA:«!A,l'.Pl,:^XLAf.PA,XCP.G) 

PLOT  nCCOUSTP.UCTlOll 
WRITn(6,55) 

URITE(G,1001*) 

IFdPlT.EO.O)  GO  TO  GSOl* 

CALL  PLOTIT(OATA2*2,MPl,l,MPl,ISKIP,t!AXPL,l,1.0> 
com  I  HUE 

UR:TE(6,2) 

Wa‘*7E(G.ia22) 

GO  TO  1*321 
com  I  HUE 


FORMAT(5IS,3F5.0) 

FORMAT! ‘1’) 

FORMAT!//, IX, 'TRUE  P.ESFOHSE  VERSES  RECOf.’STRUCTEP  RESPONSE',//) 
FORH.AT!  lUl,  SOX,  '  STARTING  SIMULATI  Ol."  ,/20X,  '  SYSTEM  ORP-R  =  ',15, 
l/,20X,*n  ♦  1  =  ',f5,///,20X,'SAMPLr!C  I’lTERVlL  =  ' ,F10.G,//,2i:<, 
2'IREM  =  ’,I5,/,20X,'ICIAS  =  ',15,/) 

FORMAT!SI5) 

FORMAT!//,5X,'i:iPUT  !♦)  AHD  OUTPUT  !•)  OF  THE  PLVIT',/) 
FORt:AT!//,3X,'  OUTPUT  <•)  AHO  RECO.HSTRUCTIO-I  !♦)  ',/) 

fop.;;at!0'jM) 

FORMAT!/////////////////////////////////////////////////////////) 
FORMAT! /,1a,  '•*••*•••••*•••*•••**»•»••*•****•*••*•**•»•**••**’'', 

1  **♦*«************«***&**  ) 

FORMAT!£F10.0) 

STOP 

OEFIIIITIOH  OF  PARVtETERS  USZP  IH  THE  SIMULATIO.H  OF  A 
LIIIEAR  OYHAMIC  SYSTEM 

X  IS  THE  CORRUPTEO  OUTPUT  SEQUEHCC 
V  IS  THE  CORRUPTEO  IHPUT  SCQUEHCE 
GAI-yiA  IS  THE  COEFFICIEHT  VECTOR 

MAX  =  ACTUAL  niMEKSlOH  SIZE  OF  2-DIM  ARRAYS  i:-  THE  DI".E!I5IO-.' 
STATEMEHT 

11  =  ORDER  OF  SYSTEM 

THE  MAXIMUM  VALUE  OF  !J  IS  !L\X/Z-1 

HPl  =  M+1,  THE  TOTAL  HUMGER  OF  5Ar..->LEn  POIIITS  IH  EACH  SEOU'-’ICE 

RUO  =  EXPECTATIOH!  W(K)*Q(K)  ) 

DELTA  IS  THE  SAMPLIHS  IHTERVAL 

IGRM  =  1  GR.AMI  IS  PERFORMED 

IPLT=3  HO  PLOTS 

IPLT=1  PLOTS  DULY  WITH  PRIHTER 

ir.lAS  =0  110  BIAS  IS  ASSUMED  PRESEUT  OH  lUPUT-OUTPUT  DATA 
IBIAS  =1  SMALL  VALUES  OF  IHPUT-OUTPUT  SI -AS  ARE  ASSUMED  PRES'.HT 
OH  THE  DATA. 


i  m 


I  s 

il 

m 

tm 

■  -91 

>  m 

•  ^ 

*  ^ 
f  ^ 


t  ^ 


END 
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SUBROUT  I  ME  CRAM » {  X,  V,  !1P1 ,  M,  DELTA, 0.3AV,  KO?T,  CAIU'A,  XLA^■.n^,  C.  Z,  AWX  , 
II  REM) 

THIS  SUDROUTIIIE  PERFORMS  TIIF  GRAKI  TECIIIJIQUE 

DI:;E!!3I0:1  X(l),V{l),G(HAX,l),ZJ;tAX,l>,GAH;.!A(l),XLVif'A(l),a(ZD), 
1DFL(23) 

DIMERS  I  Oil  GAtKZS) 

DOUBLE  I'RECISICH  GAM 

DOUBLE  PRECISIOIl  G,2.nAKItA,XLAMPA,0ELTA,nEL,PR0P,Q,QGA7 
REALMS  VARO,VARU,FAC 
REAL'S  SCALn(20),SCAL 

com.';o:j  /gkro/igkr 
coiiMOM  /iiui'ES/:i:i 
COI'-MOU  /BIAS/ 1 31  AS 

HRITE(C,10Q3) 

FOr.:nT(lHl,23X,'THE  CRAM  I  TECinllQUE*) 

JOPT  =  0  IF  DIRECT  TRAIISMISSiOM  IS  ASSUHMED 
JOPT=0 
IDLY=IREH 

iF(iREi'..:jn.o)JOPT=a 

DEL  IS  '■  !l£  HUHERATOR  OF  THE  CMOlirl  FIRST  ORDER  DIGITAL  FILTERS 
D019l=l,!l 
D£L(I)=1.0D03 
Q(I)=QSAV 
KRITE(r.,2020) 

FORMAT ( SOX, 'Q  PARAMETERS') 

CALL  PRVECCQ,:!) 

NPl=:i*l 

IIP2=:i«2 

hp:ipi«!i*i(*i 

!IPflP2=ll'H'2 

HR=!!P1-IREM 

HP1PIR=:JP1«I3£M 

DO  12  1=1, MAX 

DO  12  J=1,IL\X 

Z(l,J)=0.0D00 

VARQ=3.0 

VAR’./=0.0 

DO  500  1=1, HP! 

VARW=VARU*V(I)'V(I) 

VARa=VARQ*X{l)=X(l) 
yARQ=DSQRT  (  VARO/l’.?l ) 

VARH=nSQRT(  VARH/i:Pl ) 


CALCULATIMG  THE  G  MATRIX 
IFdCIAS.EQ.OlGO  TO  11 

iip:ip2=:i':i-»2+i 

HR=IIP1-IR£K*1 

COIlTIliUE 

D010l=l,MPriP2 

GAM(l}=C.O 

GAt^!A(l}=0.0D00 

0010J=I,IIP!IP2 

G{l,0)=0.0000 

GAK(1}=1.0 

0050K=l,iiPl 

IF(K-IDLY)2S,25,2Ii 

GAiaiACJPZ)  =0,0000 

GO  TO  2G 

FAC=1.0 

GAIIMA(HP2)=V(K-I0LY)/FAC 

FAC»1.0 
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I'S 


GAIlM^(l)  =  X{K)/t••AC 

COIITIMUE 

DOi01=X,N 

fiAlU  1  +  1)  =riAiU  t  + 1)  *a(  I ) +GAH(  I)  *Dr.L  ( I  ) 

GAMMA(  I  +1 )  ^GAMlIAt  I )  *rvLL<  I  )+GAtHlA(  l  +  l)»OJ  I  ) 

GAMHAC  1+!)I’2)=GAMHA(  I+!IP1)*[)CL(  I  )+GAt:i;A(  I♦:IP2)*0.(  I ) 

I  F  (  I  B 1  AS .  EQ.  1  )GAMItA(riPHP2  )=GAM(  I  +  1 ) 

DO  40  I=1,I1PMP2 

DO  40  J  =  I,DI’IIP2 

GCI  ,J)=G(  I ,  J)+GAI^f!A(  l)»GA(i:‘.A(J) 

COHTlIlUr. 

PRC  FORI  '  SCALIMG  Oil  G-hATlUX 
DO  735  I»1,I!P:!P2 
SCALr(l)=n$QRT(G(l,l)) 

DO  736  l=l,IIPI)P2 
DO  736  J  =  l,il,>tiP2 

G{|.J)=G(l,J)/(SCAtE(l  )*SCALr.(J)) 

WRI fE(6,1000) 

FORIiATdOX,  ---G  MATRIX—*) 

DO  55  I=1,I1PI!P2 
vmiTE(6,3)(G(J,l).J=l,l) 

FORIlAT(lX,10ni3.5) 

CONTI  HUE 
00601=2, IIPi)P2 
K=l-1 
00600  =  1, K 
G(l,0)=G(J,t) 

IF(JOPT)70,90,70 
DOGOO  =  1,IIPMP2 
00301=1, HR 

G(IIP1+1,0)=G(HP1PIRH,J) 

COHTIliUE 

HPHP2=!IPIIP2-IREI1 
00350=1, HPHP2 

'',A!.r.<HPl+0)  =  SCALtvH.nPlR+0) 

1)085  1=1, HR 

G(0,HP1+')=G(0,NP1PIR+I) 

COHTIliUE 

COHTIHUE 

CALC  GURDCT ( G ,  Z ,  XLAMO A, H, IIPtlP2 , MAX ) 

DE-SCAlC  SYHTHCTIC  COCFFICICHT  VECTOR,  XIAMDA 

DO  741  1=1,!!PHP2 

XLAIIDAC 1  )=XLA:\PA(  1  )/SCALE(  I  ) 

IF(IDIAS,EQ.J)GO  TO  43 

HPNP2=HPHP2-1 

HR=HR-1 

COilTIH',-: 

XMEAH=«UAMDA(HPHP2+1) 

IF(OOPT)120,13O,120 
HPHP2=:iPHP2+lREM 
001221=1, HR 

XLAMDA(  HPHP2-  I  +1  )=XLAI10A(HP2+HR- 1 ) 

001231=1, IREK 

XUAMDA(HP1+I)=0.0D00 

COHTIHUE 

FAC=1,0 

DO3011=HP2,IIPIIP2 

XLAMDA(I)=XLAMDA(1)*FAC 

WRITC(6,1001) 

FORI!AT(10X,'THF.  SYNTHETIC  COEFFICIENT  VECTOR,  XLAMDA,  IS') 
CAUL  PRVCC{XLAnDA,HPHP2) 

DO  150  l=l,MPHP2 
GAMMA(l)=XLAHnA(l) 

GEHERATINR  GAMMA  FROM  XLAMOA 


i  !i 


11 


'  ‘K’t-  ?  V?**  * 
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CALL  RUli.  .  .  ^EL,!l,f!AX) 

DOlbOl^l.lit.' 

GAMHA( I ) =0.0000 
1)0100>)  =  1,;!P!II’2 

GAMfiAC  I  )=GAtli;A(  I  )  +  G(  I ,  J)*XLAMDA(  J) 

COtlTKlUC 

D0200I=2,IIPNP2 

GAMI!A(  I  )  =OAM,‘iA(  I  )/GAfMA(  1) 

GAMflA(l)  =1.0000 

IF(  IDLY.CQ.OIGO  TO  172 

inLYl“inLY+l 

DO  170  I  l  =  ll)LYl,»Pl 

1=HPHP2+1-I  I 

GAIit:A(l  +  IDLY)=GAHMA(l) 

no  172  l“l,inLY 

CAIiMA(  I +m>l)  =0.0003 

COUTIIIUE 

CALCULATK^G  Tllf  EQUIVALEtlT  CCriTIIIOOUS  OESCRI PTIOM 
CALL  l7.T0S(GA!:.>lA,:l,nCLTA,K0PT) 
imiTE(G,1033) 

FURIMT(/// JX,  100(111- ),/,  IX,  100(i:i-)) 

rctur:i 

END 


SUiiROUTlifE  GKRnCT(X,Y,XLAMnA,!ltl,:i,l'AX) 

REaL*8  X(ltAX,l),Y(MAX,l),A,5,C,n,C,DCT,CCC(20,20),XL\ttPA(l) 
IDTEGER  !n:!l(2,20) 

COlillOli  /OKRO/IGKR 

IGKR=0  USE  IS  MADE  OF  TIIF  FIRST  ROW  OF  ADJOINT 

1  DIAGO-IALCIEGATIVE  ENTRIES  SET  TO  ZERO) 

2  ABSOLUTE  VALUE  OF  DIAGOtlAL 
DO  0  l=l,N 

DO  G  J=l,ll 
Y(J,I)“X(J,I) 

A=1.0D0 

DO  45  1=1, Jl 

B=0.0D0 

L=l 

tl=l 

FIND  LARGEST  ENTRY  A(L,tl)  IN  LOWER  DIAGONAL  SURHATRIX 

DO  13  J=I,N 
DO  13  K=I,U 

IF(OADS(Y(K,d)).LE.B)GO  TO  13 
B=DA3S(Y(K,J)) 

L=K 

H=J 

CONTINUE 


INTERCHANGE  ROWS 

IF(L.EQ.I)G0  to  24 
DO  25  J=1,N 
C=Y(L,J) 
Y(L,J)=Y(I,J) 


23  Y(I,J)»C 

C 

C  IMTCRCIIAdCiC  COLUMNS 

C 

21)  ircn.ra. I )c«o  to  23 

no  28  J“1,U 

Y(J,M)»Y(0,l) 

28  Y(J,I)»C 

C 

c  OEGtit  sucr.p  colunus  to  tuc  right 

C  ARRAYS  tlUI-Ul,.)  ,mitU2,.>  KELP  RCCORO 

C  OP  RO\l  AND  COLUl"!  I  SUr.RCUANOr.S 

C 

23  IIUM(1,I>»L 

HUM(2,1)-M 
B*Y(I,I) 

Y{M)«A 
no  1(2  J«l,!l 
lP(J.r.a.  t  )GO  TO  ((2 
C— Y(I,J) 

Y(i,o)“0.ono 
no  1(1  K«i,ti 

t)»Y(K,l)*C 

E-Y(K,d)«r.*n 

c  tr(oAiiS(r.).LT,i.on-io*nARS<o))n"0,ono 

41  Y<K,J)<-t;/A 

42  CONTI  IlUC 

43  A-O 
C 

C  RCSTORC  COLUMNS 

C 

no  SS  1-2, N 

d-lul-1 

K«NUM(2,J) 

IF(K.CP,.J)00  TO  52 
DO  51  L»l,ll 
C-Y{K.L) 

Y(K,L>«>Y(J,L) 

51  Y<vl,L)»C 

52  K-NUM{l,d) 

C 

C  RESTORE  ROWS 

C 

IP(K.CQ.O)GO  TO  53 
DO  57  L«1,N 
C»Y(L,K) 

Y(L,K)-Y(L,J) 

57  Y{L,>J)»C 

58  CONTINUE 
OET-A 

C  •***«**•  SET  IPRIIIT  ***** 

IPRINT-O 

IPRKIT"! 

IFdPRINT.NE.DGO  TO  111 
WRITE(0,101) 

101  FOmiATC/lX.’DET  01  GRAM  MATRIX  IS  ’) 

CALL  PRVtC(OET,l) 

WRITt;(r),102) 

102  F0RMAT(/1X, 'AOdOINT  MATRIX  IS') 

CALL  PRMAT(Y,N,N,MAX) 

001001=1, N 

001000  =  1,11 

CCC{|,0)»0.0000 

D0100K»1,N 

100  CCC{ I ,0)=CCC{ 1 ,0)+X( I ,K)»Y(K,0) 
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I 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


JO 

950 

993 


1000 

1001 


200 

C 

C 

C 

C 

C 

C 

c 

c 


15 

7 

UG9 
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SUniiOUTItfi;  ZT0S(ri,A/i,0F.tTA,l7.1S) 

coHHOH  /uunr.u/:ni 

COHVrUSlOU  OP  A  niSCRCTE  TIIIC  SYSTEM  IKZ)  TO  A  CONTIIIUOUS  TltlE  SYSTEM  ll(S) 


1I(Z)»(A(1)  ♦  A(2)*zrTA  « - )/(l  ♦  n(2)«ZETA  ♦....) 

ZETA  »  1/Z 


ius)  =  (A(i)  ♦  A(2)‘r.  ♦  .  ♦  A(:i*i)*s**:i)/nE:iot‘. 

DE!)0M»IU1)  ♦  G{?)*S  ♦  . ♦  D(!l+1  )*G**I1 

R(l)  »  1  ALWAYS 


OlllEIISIOIl  R(9),A{0),Tr-.P(20),n'U20),lU(20),C!l(?9),C\(?1),CA\(?.0), 

^COIU'LExilG^CAlcAAicAi;cr.Icfi'c6TlIcW2,COMT,r\r,Al,A2,Rl,02,A\l,IO-^ 


LU,  i  ,  v.r  « 

REAL*.".  H,A,TEnP,RR,RI,nELTA 

C0IIT  =  0. 01100 

lORP’IZTS 

NPl*:)*! 

IIUP1<-.|!;)<  1 
Al"0.0!ir. 

U1"0,0!)0 

DO  JO  1  =  1.  IIP! 

ir(l  .LE.:!I1PI)AI“A1*A(  I) 

ni=r.i*r.(  i ) 

WRITE(u,039)  RR 
rOP.!!AT(  lux, ':»!  =  ’,15) 
rORMAI (///) 

WRITE <0.999) 

WHITECG.IODj) 

FORim  ( '  7.-llO!'.AI  IJ  OroOtimATOR' ) 
CALL  i’.{VEC(!:,::ri) 

WRITECo.lOOl) 

rORMAU'  Z-l)0:iAI!l  tUM'.ERATOR') 
CALL  PRVrC(\,:iPl) 

IFdZTS.rO.O)  GO  TO  919 
I  F(  I  ZTS.r.0. 1 )  GO  TO  230 

irciZTs.Ea.z)  go  lo  zso 

ir( IZIS.EQ.J)  GO  TO  200 
irdZTs.rQ.-'i)  go  to  zso 

CORTII.'UE 


LOGAR  I  Tlllll  c  TRAIISEORI’.ATI OR 


WORK  0:i  UUlirRATOR 


IP<NN.r.!K0)G()  TO  4G9 

CALL  POI  lUC  •„TCnP,ri:i.RR,RI,IER) 

no  15  l  =  ',l!ll 

CAd  )  “tici '.I’EX  ( RR d  ) .  R I  d  ) ) 

DO  7  i"i,n:i  , 

CA(  I  )“(♦!  .0/01  LlA)*CnLOG(CAd)) 
irdllM  (l.!ll  0010471 

com  I  laiE 

DO  47.1  l=:i!IPl.MPl 
CAA(  1  )-0.0!'3 
CAd)“3.0P'J 
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47 1  CONTIHUC 

IF(l<H.Ca.O)CAA(l)»l.onO 

C 

C 

C  NOW  THE  FIRST  IIN  EHTRIES  OF  CA  COKTAIM  THE  S-nOMAIII  ZEROES  OF  NUMERATOR 
C  Aim  THE  REAHIHING  ENTRIES  ARE  ZEROED  OUT. 

C 

IF{NH.HE.O)CALL  POLCOfHCA,CAA, 0,H) 

C 

C 

C  S.'ORK  ON  DENOMINATOR 

C 

C 

919  CALL  POLRT(B,TEI',P,!l,RR,RI,  ir.R) 

OOIG  1-1, N 

cR(i)«nci;pi.x(RR(i),rvi(i)) 

1C  CF(l)=1.0mO/OR(l) 

909  '.IRIIEIC.IOOZ) 

CALL  PRC\'EC(Cr,ll) 

IFdZTS.EQ.O)  GO  TO  900 
255  OOC  I -1,11 

0  CRt  I )»(-!. 0/r.rLTA)«CnLOG(CR(l )) 

WRITC(f.,2!iO) 

'240  rORMAT<'  lOGARITHM.tC  TRVISFORHATIO'l* ) 

WKITF(G,!>99) 

WHITE(C,2000) 

2000  rORIIATC  POLES  IN  S  OOIIAIN') 

CALL  PRCVCC(CR,M) 

0050001-1, M 
5000  CR(I)-'CR(I) 

CALL  POLC0N(CR,Cn,0,N) 

C 

C 

C  ADJUST  DC  GAIN  CONSTANT 

C 

C 

A2=CAA(1) 

B2»CB<1) 

rAC=<Al/i>l)dl!2/A2) 

DO  e05  I-1,NNP1 
005  CAA(l)«CAA(l)«rAC 

GO  TO  2010 
C 
C 

C  PELAYFP  PULSE  INVARIANT  TRANSFORMATION 

C 

C 

C  SHIFTS  NUMERATOR  COEFFICIENTS  FOR  DELAY 

C 

C 

250  CONT-A(l) 

DO  300  I -I,!! 

303  A(l)»A(l*l)-r.ONT*n(l*l) 

A(IH>1)  =  0.0 

400  CALL  P0LI:T(:'.,TC.’!P,H,RR,RI,  lER) 

DOoll“l,!J 

CI!<l)=i!CMPLX(RR(l),fil(l)) 

Cl  cr(i)=i.n:>03/cu(t) 

WRITE (0,1002) 

1002  rORHATdX, 'THE  POLES  OF  THE  Z-DOMMN’) 

CALL  PKCVEC(CF,N) 
t 
C 

C  PARTIAL  FRACTION  EXPANSION 

C 

C 
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m>^-m 

flaatljS- 


D03l=l,» 

CO!11=1.0POO 

COH2=O.ODOO 

D0I*J“1,N 

C0M2  =  C0!I2*CR{  I  )  +  A(N-J+l) 
IF(I-J)5,I(,5 

COII1=COII1*(1.0000-CR(I)«CF{J)) 

COHTIJIUE 

CA(I)=C0!)2/C0II1 


TRAHSFORMATIO:i  OF  DCIIOMI JIATOR  AtIO  NUMERATOR 


D021=l,:l 

CR(l)=CDIOfi(CR(l))/!!ELTA 

CA( I ) =CA( t ) ‘CRC I ) / ( 1. ODOO-CF( I ) ) 

CONTI  HUE 
cA(iiPi)=o.onoo 

IIRITr(G,2Ul) 

rOUMATC  DELAYED  PULSE  TRAMSFORHATf  OH’ ) 

»RITn(G,999) 

l)RITG(G,109!i) 

FORMAT  (■  NEGATIVE  OF  THE  POLES  IM  THE  S-DOnAi:i') 
CALL  PRCVEC(CR,N) 

HRITE(G,100J) 

rORMATdX/ NUMERATOR  CONSTANTS  OF  FACT ORl ZED  IKS)') 
CALL  PRCVEC(CA,::) 

CALL  POLCON(CR,CB,0,N) 

D071I=1,NP1 
CAA(  I  )  =  '.). ODOO 
009K=1,II 

CALL  POLCON(Cil,CFl,K,N) 

D09J=1,I1 

CAA(J)=CAA(J)+CF1(J)*CA(K) 

CAA(NPi)=y.onoo 

CONTINUE 

DO  1)59  1=1, NPl 

CAA( I )=CAA( I )+CONT«CB(  I  ) 


WRITING, 1005) 

FORMAT! '  S-OOMAIN  DENOMINATOR') 
CALL  PRCVEC(Cr.,NPl) 
HRITC(6,100G) 

FORMAT!'  S-DOMAin  NUMERATOR') 
CALL  PRCVCC!CAA,NP1) 
no20i=i,ripi 
B!l)=cr.!l) 

A!I)=CAA!I) 

RETURN 

END 


"'f  ■ 
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SUnP.OUTI!IC  POLRT{XCOF,COF,;'.,llOnTR,ROOTI,  l-n) 


COMPUTPS  TUP  nPAl  AND  COMPLF.X  ROOTS  OF  A  RPAL  POLYROf.l  \1. 

ncscniPTiO!!  or  parai!Ctfrs 

xcoF  -vncTOPv  or  corrricirriTs  of  Ti:r  polyoo."! m. 

ORPKRPP  FROt;  Sf.AI.LF.ST  TO  I.ARGCSI  PO-TR 
COF  -\J0RKi:!0  vnCTOR  or  LFTCTII  tl+1 
M  -OROPR  OF  POLY.TOtllAL 

ROOTR-RPSULTATT  VCCTOR  OF  tPrjf.TII  fl  OOriTA I ‘1 1 ‘Jr.  RPAL  ROOTS 
OF  TMP  POPYTOlilAI. 

ROOTI-RPSU!.TA!!T  VPCTOU  OF  LPURTH  tJ  COtITA I  : II.T.  TUP 

cop>i!P5PO!ini::5  ir-AoirmY  roots  of  tmp  poi.y’io;-i \i 
lER  -pi;i:o;:  coop  i.tierp 
tPPv-O  fio  ERROR 
IPRr.l  (1  LCSS  TIIA:1  OTP 
IPR=2  M  GRPATPR  TMA’I  JG 

lEP--^  UTAIllE  TO  I>PTPP.!!I!1P  ROOT  WITH  SOO  I TTPRAT I O.TS 
Oil  5  START  I  Ilf.  VALUES 
tPR.-l|  moil  ORPPR  COEFriCIPIIT  IS  ZERO 

DII-.PIJSIOII  XCOF(l),COr(l),ROOTR(l),ROOTI(l) 

nOur.LP  FRECISIO:;  XO,YO,X,Y,XI>.R,YPR,OX,UY,V,Yl.XT,U,XT2,YT2,StlMSn, 

1  |1X,DY,TPI!P,  ALPHA,  XCOr,  cor,  ROOTU,  ROOT  I,  PR1,CR2,XSS,X5,YS5,YS,T0L 

Lir.lTPn  TO  3GT!I  ORDER  POLYIIOMIAL  OR  LCSS. 

FLOATKIG  POIIIT  OVERFLOW  HAY  OCCUR  FOR  HIGH  ORDER 
POLYIIOIIIAI.S  RUT  WILL  HOT  AFFECT  THE  ACCURACY  OF  THE  RESULTS. 

METHOD 

IIPWTON-RAPIISOT  ITERATIVE  TCCIIIIIOUF..  THE  PITAL  ITPRATIOHS 
OH  EACH  ROOT  ARC  PERFORIIPO  USItlO  THE  ORIGITAL  POLYTOf.l A!. 
RATHER  THAU  THE  REDUCED  POLYHOHIAL  TO  AVOID  ACCUr.ULUT) 
ERRORS  III  THE  REDUCED  POLYHOHIAL. 

ER2"l,0n+5O 
TOL^l.OO-S 
IFIT=0 
II =M 
IER=0 

tF(XC0F(!l+l))10,25,10 
IF(H)  15,15,32 

SET  ERROR  CODE  TO  1 

ICR=1 

IF(IER)200,201,200 
WRITn(o,203) IPR 

FORI!AT(  IX, 'ERROR  CALLED  FROM  POLRT,  lER  »  ’,13) 

RCTURH 

SET  ERROR  CODE  TO  4 

lER-H 
GO  TO  20 

SET  ERROR  CODE  TO  2 


lER’a 
GO  TO  20 

IF(II-3G)  35,35,30 

HX-II 

IIXX«=H*1 

N2-1 

KJl  -  ll+l 
DO  40  L"1,KJ1 


5A 


o  o  o  r>  o  o 


tiO  C0F(I1T)=XC0F{L) 

SET  INITIAL  VALUES 

US  X0=. 00500101 
Y0=0. 01000101 

ZERO  INITIAL  VALUE  COUNTER 

IN=0 
50  X-XO 

INCREMEMT  INITIAL  VALUES  AND  COUNTER 

X0=-10.0‘YO 
Y0=-10,0*X 

SET  X  AND  Y  TO  CURRENT  VALUE 

X«=XO 
Y=YO 
IN=lii*l 
GO  TO  59 
55  IFIT-l 
XPR=X 
YPR”Y 
C 

C  EVALUATE  POLYNOMIAL  AND  DERIVATIVES 

C 

59  ICT=0 
GO  UX=0.0 
UY-0.0 
V  «0,0 
YT=0.0 
XT=1.0 
U«C0F(IH1) 

IF(U)  G5,130,G5 
G5  DO  70  l»l,N 
L  »N-I+1 
TEMPaCOFd.) 

XT2»X*XT-Y«YT 
YT2=X*YT+Y*XT 
U»U+TEMP«XT2 
V»V-»TCMP*YT2 
Fl  =  l 

UX«'UX*FI*XT*TEi;P 
UY“UY-FI*YT«TCMP 
XT=XT2 
70  YT»YT2 

SUMSO.=UX*UX+UY*UY 
IF(SUMSQ)  75,110,75 
75  OX=(V*UY-U*UX)/SUMSQ 
X=X+DX 

DY— (U*UY+V*UX)/SUIlSq 

Y=Y+0Y 

XSS=X 

YSS=Y 

IF{YSS.Ea.0.0D0)YSS»1.0O0 

IFCXSS.EO.O.ODOlXSS-l.ODO 

ERl-OAr>S(DX/XSS)*DAnS(DY/YSS) 

IF(ER1.GT.ER2)G0  TO  73 

ER2=ER1 

XS=XSS 

YS=»YSS 

73  IF(ER1-TOL)100,30,30 
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c 

c  STEP  ITER.UIOIl  COUNTER 

C 

80  lCT»tCT+l 

IF(lCT-500)  C0,S5,S5 
85  IF(IFtT)100,90.100 
90  IF(ll!-5)  50,95,95 
C 

C  SET  ERROR  CODE  TO  J 

C 

95  IER-3 
X«XS 
Y"YS 
ER1“EK2 

100  no  105  L=i,:ixx 

ht«i:ji-L'*i 

TEUP'tXCCFd'.T) 

xcor(t:T)»coF(L) 

105  COF(L)»TEf.P 
iTCi'.p=:i 
N=tJX 

!!X=ITr.UP 

IF(lFn)  120,55,120 
110  IFdFlT)  115,50,115 
115  X^XPR 
Y-YPi; 

120  I  FIT-0 

122  lF(nA:'.S(Y)-1.0n-S*nAr.S(X)  1155, 125, 125 

125  All'li.\-X«X 

SUrSO=X*X*Y*Y 
i)»II-2 
GO  TO  11(0 
150  X-1.0 

iix«:i,\-i 

tixx«:!xx-) 

155  Y-Ci.O 
sursn-’C-.o 
AiPi:\«-x 
1 

no  COr(?)<^CCF(2)'*AEP:i\*COF{l) 

lUS  no  150  L“2,'J 

150  CO."(L*l)=COF(L+l)  +  Al.P!IA*COF(l.)-SU:iSn«COr(L-l) 
155  ROOTI(;!2)=Y 
Rom!:{i!2)=>: 

lFCr.ni.GT.TOt)\ir>ITr(G,55ii>:J2,ERl 
55!i  ronn\T(ix,'CRnon  o:i  ',15,’tm  root  is  *,010,5) 
Ln2-'i.on«:.o 
!!2=:!2*1 

IF(SUI!SQ)  103,105,150 
100  Y'-Y 

sur.sa^o.o 
GO  TO  155 
105  IFCI)  20,20,45 
E!in 
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SUnROUTII.'E  PLOTIT(X,tlF,llPT,ll,l2,IJ,i;\XX,ISC,SC,\Ln) 

DllirtlSIO!!  X(ll.\XX,MF),irL(07>,IP(07),l(;H(9),IS(n),X!'X(0),XI!'!(9),F(9) 
1) 

RrAL«3  X.”\X,Xt'l!J 

DATA  IC!l/lH«,lH*.\HJ,ltii:,lllS,lHG,lH7,l!l5,lll9/ 

DATA  llll,>;i,lliG.MlL/lHt,III  ,l(l>,l:K/ 

D0201l»l,97 
201  IPCD'III 

002001-1,8 
0-12*1+1 
200  IP(J)-II!I 
IP(1)«I!!I 
IF(liF-J)15. 15,17 
15  NP-lir 

GO  TO  IS 


!  i 
!  ^ 


'  ^ 


17 

IS 


5 

1 

111 


117 

US 

116 


115 

111) 


50S 

507 

5 


51S 


!1P=5 

ooiiij»i,::f 

x:ii:;-x(i,j) 

XHAX-X(1,J) 

DOii“2,!:rT 
IF(Xf.\X-X(  l,0>)'',2,2 
iF(xmii-x(i,j))i,i,u 
XI!1!!-X(I,J) 

GO  TO  1 
X!lAX«X(t,0) 

coi.'Tni'i: 

XfX(0)-X!-.\X 

xti!J(0)-x:'i'! 

UFF-HF 

iFdsc.-JiM.OR.x'F.n^.ur.o  to  iic 

xi'.AX-xi;x(i) 

xmn»xi'.N(i) 

001170-2, ’jr 

IF(X!-AX,UT.XIiX(0))X!'AX-XnX{0) 

I  r  ( xn  I ,  r.T .  xiri  ( 0 ) )  x!’  I  vm:.'  ( o ) 
cn-n  inur 

001130  =  1, 'IF 
XI!X(0)=XI11X 

xi- HUoi-xm-j 
NFF-l 

n  (scALr.!!M.i)r.o  to  m 

001150-1, :irr 

xx«o.5«(\i.x(.i)  +  (i  +  scAL!:)+x".:uo)*(i-scur>) 
xivi(o)=o,  ,  +  (x:-.x(o)*(i.o-scALC)+xp:uj)  +  (i.9+r.CM ru 

XI!X(0)=XX 
0011150=1, X'pr 

xiiAx-xr.x(j) 

xi',i:i»xi::J(o) 

xx»oAr.S(:;i;AX) 

xii- OABS(::f.i:i) 

IF(XX.LT,x:!)GO  TO  505 
IlAX-1 

XSAV-XX 
GO  TO  507 
}1AX  =  0 
XSAV=X:i 

1  r  ( X! :  \x-  x:  I  :j-  1 .  on-  G  *  xs  AV )  5 , 5 ,  G 
XPAX  =  XI’.AX+XX+1.0n00 
XMI!J=Xlllil-x:i-1.0n00 
xx=DAr.s(xi;AX) 
x:)»OA:'.s(xm:i) 

IF(XX.L1 .XUIGO  TO  513 

JlAX-1 

XSAV-XX 

GO  10  G 

!lAX-0 
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XSAVXK 

,  X«=At.Or.lO(XSAV) 

JJ=X!I 

XX=JJ 

lF{X!J-XX,!IK.0.O.ANn.XSAV.LT,1.0)JJ“JJ-l 

Tnw=io.ov»jj 

K=XSAV/(Tnil*0.0) 

ird'.AX.CQ.DT.O  TO  503 

IF(Xm:l.l.T.0.0)t>-(K*l) 

XX»XMAX,TCil 

KI!I!I>>!>3 

X:i=Kf.l!! 

505  XiJ'X.’l+O.O 

tF(X:J.LT.XX)r.O  to  505 
W*AX“Xtl 

I  F  ( KCAX*  Ktl  I  :i .  tT.  0 )  K(iAX=-Kltl N 
GO  TO  1112 

503  IF(Xi:AX.GT.0.0)K=-(K*1) 

xx«xni»/T£:ii 

KIIAX— K*S 
XIJ=K‘*.AX 

50!)  XIJ=-Xtl-3.0 

if(xii.gt,xx)go-to  so; 

KHIN'Xtl 

I  F  ( KHAX*  K!'.  I  iMT .  0  )  i:iJ  1 :1=-KI1AX 

1112  Xt'.X(J)=KI.AX*TE:! 

1113  X!lN(J)“i;!n!l«TEii  ^ 

IF(  ISC,:!C.1.0R.HF.EQ.1)G0  TO  119 

n01200=-2,«F 

XtiX(J)<=XltX(l) 

120  Xt!!J(J)«X!:!l(l) 

119  D01120=--l/<Fr 

r(i)=x:;!UJ) 

XIII II  >=(X!;X(0)-XM!»(J))/8. 0030 

0071=1,0 

7  r(i+i)=F(i)+xin!i 

IF(F(l)*r<9).tT.0,0)F(5)»0.0 
\miTF.(G, 100)0, F 

r0!ir.AT(5X,  '  ruilCTION'  ,  I  2,7X,3E12.I)) 
l.'RITC(6,106)ICII(0),IP,II’ 

FORMAT ( 5X,  '  SYliUOL  '  ,  IX,  A1,1SX,97A1,/ ,33X,  97A1 ) 
0051=1,97 
IPL(l)=in<!) 

00113I=1,?JF 

F(  I  )=9G.0030/(XHX(  I  )-XllH(  D) 

OOGOOI=1,IIF 
IS(n=l 

IFCHP.EQ.DASSICH  1001  TO  HPPP 
IF(IIP.CQ.2)ASSIGII  1002  TO  UPP? 
IF(IIP.nQ.3)ASSlGK  1003  TO  HPPP 
0093=11,12,15 
KI1AX=1 
KIUK=97 
00101  =l,:iF 

K-F(n*C-(3,n-XI«l(l))*1.5 
IS1=IS(I) 
ir(K)662,GG2,o621 
6G21  IF(97-K)G331,G53,G33 

C331  K»97 

ICCC»niG 
GO  TO  Gl)U 
G62  K»1 

ICCC=1!IL 
GO  TO  CI4I1 
635  ICCC»JCII(t) 
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112 

lOG 


113 

GOO 


I 

i 

f 


f 


s 


IF(ISt-K)G011,C01,G01 

D0G02L=JSI,K 

IPL(L)“ICCC 

iF(KMi!i.riT.isi)Kr.i:)=isi 

IF(KtlAX.LT.K)KMAX=K 

GO  TO  u03 

D0C3ItL=K,  ISI 

l?L(L)=ICCC 

IF(Km:i.GT.K)l'tlltJ=K 

I F ( KMAX .  LT .  I S I  )  K'!AX=  I S I 

ISC  I )=K 

CONTINUE 

GO  TO  lIPPP, (1001, 1002,1003) 
IJRITE{C,102)J,X{J,1),IPL 
FOmUTClX,  |I|,E1?.!»,1GX,97A1) 

GO  TO  IS 

WRITE(G,103)J,X(J,1),X(J,2),IPL 

FOR!1AT(lX,l4,2E12.1i,l|X,97Al) 

GO  TO  IS 

WRITE(G,104)J,X{J,1),X(J,2),X(J,3),IPL 

FORMAT{lX,lli,3E9.3,lX,97Al) 

0019l=Km!l,KllAX 

IPL(I)=IP(I) 

CONTINUE 

WRI7CfG,105) 

FORImATCIUI) 

RETURN 

END 


SUBROUT  I  HE  RESPON  (  X,  V,  N,  GANIIA,  XLAIlDA,rPl ) 

DIMENSION  X ( 1 ) , V ( 1 ) , GAHMAC 1 ) , XLAMPAC 1 ) 

REAL*S  XSAV,GAF'J1A,XLAMDA 

HH1=N-1 

HP1=N+1 

HP!IP1=N*N*1 

HPNP2=:i+N+2 

DO  19  i=i,np:ipi 

XLAMnA(l)=0,0n00 

XSAV=0.0P00 

DO  20  K=l,nPl 

IFCl.EQ.DGO  TO  25 

BO  21  1=1, NIU 

J=NP1-I 

X  L  AMDA  (  d  )  =X  LAI  lOA  (  J- 1 ) 

CONTINUE 
DO  22  1=1, H 
J=HPHP2-I 

XLAMOA(J)=XL.VinA(J-l) 

XLAHDA(1)=XSAV 

XLAi5nA(nPl)=V(K) 

XSAV=O.ODOO 

DO  23  i=i,np:ipi 

XSAV=XSAV-GArJW(  I  ♦1)*XLAI1DA(  I  ) 

IF(0ADS(XSAV).GE.l,0O10)XSAV=0.bD00 

X(K)=XSAV 

RETURN 

END 


SOBROUT I riE  ERROR( XREC, V, GAra W.HPl , K, XLAMDA, XORG ) 

DIMEIISIO.’I  XREC(l).V(l),XOKG(l) 

DltlEI.'SIOU  VVV(20) 

REAL«3  GAiaUCD.XLACDACl  ).AVr.W,SUflV2.AVGQ 
CALL  RESPO:i ( XRCC . V, H, GAtLVA,  XLAMDA^MPl ) 

AVGU=0.0D00 
SUr.V2  =  0.0P0 
D02f.t»l,i;Pl 

SUttV2=SU!;V2+XORG(  I  )«XORC(  I  ) 

AVGQ=XORG(l)-XSEC(l) 

2  6  AVGW=AVG;;*  AVGQ*  AVGQ 

AVC>.;=AVG'..7SU!1V2 

avgq=dsqi:t{avgw) 

AVGQ=100.0»AVGQ 

AVG'.!=100.3*AVR-.; 

WR I TE ( 6 , 2  7 ) AVGW, AVCQ 

27  for!;at{1x, 'PER  cs:!T  i:ca;;  fo-.jer  error  of  reco:>structio-:*,f3.5,///, 

11X,'PER  CLIiT  OF  SQUARE  ROOT  OF  POUER  ERROR  I'l  RECOSTRUCTIOrj*  F2  3) 
RETORtl 
E!iD 


sunp.ouTir.T.  BinLnA(A,Q,t)EL/:,!iAX) 

REAL«3  A(:!AX,l),Q(l),nrL{l),PRnC 

iip;jP2=!:'*:i*2 
A(l/lPl)  =  l,3n08 
PROf)=l,0030 
D03i2::=i,:! 
i>‘:jpi-k 

PROD-PROr/DELd) 

Ad,  I  )=ri;(i:) 

312  A(K«l,:iP3)  =  1.0nO0 

005131  =  2, ::pi 

D03l3r,»=i,:! 

j=:!Pl-K 

513  A(I,J)  =  (A(  l,J+l)-Q{J).A(I-l,J*l})/t>:L(Jl 
D051lt|-=1,:'P1 
rio5i!<j=i,:iPi 
A{l,J<?!m=O.OPOO 
A{I  +  :iri,. 11  =  3. 3003 
31li  Ad«!lPl,J*,IPl)=A(l,J) 

HRI7E(G,1035) 

1005  FORI ; AT  ( 3  X ,  '  A-! AT R I X '  ) 

CALL  pr.:;AT ( a, :jp!;p2 , :jp;ip2 , pax ) 
retur:i 

E.'ID 


SUr.ROUT  I  HE  PRI'AT (  A, R, II,  !!!;AX ) 

OOIICLE  PRECISION  A 

TIMS  sur.ncuTM.'r  ouTPins  pousle  precjsioii  double  niRPRsioiirn  array 
dii;e:;sio-:  a  cum,!) 

HRITE(G,1) 

D02I=l,:i 

WRITE(5,3)(A(I,J),J=1,:;) 

FORI'.ATdX,10R15.5) 

WRITE (3,1) 

HRITECC,!) 

for!;at(/) 


ft 


1 


SUBROUTUIC  PRCVEC(A,tl) 

THIS  SUItROUTItlE  PRINTS  OUT  A  COMPLEX  SItlOLC  D IMENS lOlinn  ARRTY 
A  COMPLEX  tIUMIlER  OF  THE  FORM  A  ♦  D  J  IS  OUTPUTTEO  I'l  THE  FORM 
(  A,  R  d)  WHERE  J  •  SOUARE  ROOT  OF  -1 
DIMENSION  A(l) 

COMPLi;X*10  A 
WRITC(6,2) 

WniTn(G,l;(A(l),l«l,N) 

FORMATdX  lH(,ni7.10,lH,,ni7.10,3H  J)) 

WRITE(G,2) 

WRITE(G,2) 

FORMAT(/) 

RETURN 

END 


sunROuriNE  prvcc{a,n) 

THIS  SURROUTINE  OUTPUTS  DOUBLE  PRECISION  SIlir.LE  DlttCIlS I O'lFO  ARRAY 
DIMENSION  A(l) 

DOUBLE  PRECISION  A 
URITi;(0,31) 

WRITE(C,l)<\(n#l»l»N) 

FORMATdX, 101)13. 5) 

WRITE(C,31) 

WRITE(6,31) 

FORMaT(/) 

RETURN 

END 


SUBROUTINE  POLCOIUC, R2, K,!l) 

A  POLYINOMIAL  CONSTRUCTION  PROGRAM  NEEDED  FOR  ZTOS 

DIMENSION  Cd),R2d) 

COMPLEX*]C  C,R2,COMP 
RCAL*8  DCd) 

EOU I  VALENCE  (COMP, DC) 

NPl-IHl 

D010I»2,NP1 

R2d)»o.orno 
R2d)-1.0D00 
1)04 1  "1,11 
COMP-Cd) 

IF(I  .EQ.!:.OR,(DCd).Ca.0.0D0.ANO.nC(2).rQ.0,0O0)>r.O  TO  •, 
D02JJ"1, I 
d-^i-  )J»1 

l..'Cd*l)-R2(d*l)«Cd)  +  R2(J) 

R2d)»R2d)«Cd  ) 

CONTINUE 

RETURN 

END 
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HPMINV  PROGRAM 
LISTING 


200 


fUOUltAH  llAf.l:  IIPMIIIV 

00U8LE  I'RECISIOIl  IIATRIX  IfJVCItSIO!!  PACKAGE 
Dountr.  PRCCISIOIl  AA,A,P,Q,r.,C,D\,G 

oouian  puKClsio!)  hires, rDin^rsnip, inn 

OOUliLE  PRECISIOIl  F  \C, TRAC, FFAC, ERROR 
Oltir.ilSIOn  AA(20,20),A(20,20),P(2.),20),Q(20,20),!U2'),20), 
1C(20,20),OA<20, 20), 0.(20,20) 

COMMOU  /SOL/FRAC 
coitnotl  /PR/ I  PR  HIT 
cor.nou  /PRon/PDiF,anir 

IIAX»20 

T:mr.s»i5.onoo 

WRITC(G,100) 

REAOCS  ,2li0)!l,  I  SP,  I  SQ,  I  PERT,  I  Stin,  ICOim,  ICOP.OT,  1  PRIt'.T 
WRITE(r.,  11)0)11, 1  SP,  ISO.,  I  PERT,ISLIf),  ICOIID,  ICORCT 
FRAC'l.ODOO 

IF{IPnRT.0.T.tl.Atin.lSLin.EQ.l)FRAC=1.0D-0i) 

WRITC(G,100) 

DO  200  1"1,11 

REA0(5,250)(AA(1,J),J=1,H) 

TEH»1.0n-00*AA(l,l) 

t  F ( 1  conn . na.  1 ) AAC 1 , 1 ) =aa( i , i )+teh 

COllTIllUE 

CALL  r,E()UAT(!l,nAX,AA,A) 
lF{IPRIin.ECl.l)URITF.{G,lSO) 

I  F  (  I  PR  1 IIT .  EQ.  1  )CALL  PRl!AT(:i,llAX,  A) 

SCALIllG  I'.ETHOD 

PDIF“0. 0.000 
QDIF=0.0n00 

ROW  SCALIllG 
IFdSP. HE. 1)0.0  TO  10 
1S=1 

CALL  MSCALE  (t!,MAX, TllRES,  1 S,  PO I F, A. P) 
IF(lPRIllT.2Q.l)-.miTr.(G,127) 

I  F  ( I  PRI IIT  .  EQ.  1  )CALL  PRHATCl.MAX.A) 
ir{lPRlllT.EQ.l)\miTC(C,101) 

1 F  ( I  PR  I  NT .  EQ.  1 )  CALL  PRllAT  ( :l,inx,  P) 

CALL  OI!lV(H,nAX,P) 

IF(IPRI|IT.EQ.1)WR|TE(G,102) 

IFdPRINT.EQ.DCALL  PR'1AT(II.MAX,P) 

COllTIllUE 

COLUIlll  SCALIllG 
IFdSQ.IlE.DGO  TO  20 
tS«0 

CALL  MSCALE  (11. IIAX, HIRES,  IS, QniF,A,Q) 
IF(IPRII1T.CQ.1)WR1TE(G,103) 

1 F  (I  PR  I  IIT .  EQ.  1 )  CALL  PRIIATdl.MAX.Q) 

CALL  I)I1IV(II,MAX,Q) 

IFd  PRINT.  EQ.l)imiTE(G,10ii) 


63 


ii 


1* 

I 


lF(IPRItlT.r.a.l)CALL  PHMAT(N,tUX,Q) 

20  COHTCIUE 

C  Etin  OF  SCALKIC 

IF(IPRlNT.Ea.l)WRlTC<G,]'15) 

I  r(  IPRlilT.ra.DCALL  PRHAT(!I,HAX,A) 

CALL  MEOU.\T(:i,MAX,A,C) 

if(ipert.i;q.o)c.o  to  oo 
c 

C  PERTURIi  MATRIX  A,  C»A+(FAC)«(niAr.  A)»A*FAC*(DA> 

C 

FAC-1. 3D-0it 
00  30  IIIC-1,1 
FAC-FAC*! .10-02 
\;RITE(G.100) 

WRnC(G,G00)FAC 
00  05  1-1,11 
00  35  J-l,:i 
0A( t,J)-0.0000 
B{I,J)-A( l,J) 

IF(  IPERT.GT.IDGO  TO  S2 
IFCI.IIC.IPCRTIGO  TO  SJ 

82  COliri.TUL 

IFd  ,Ea.J)a(l,d)-A(l,J)  +  FAC*A(  .,o) 

IF(I.Ea,d)OA(l,J)«A(l,d) 

83  COIITIIJUE 
C(l,J)=n(l,J) 

85  COIITIIIUE 

IF(IPRIIIT.EQ.l)  »RITC(G,130) 

IFdl’Rltn  .no.DCALL  PRI1AT(1I,MAX,Cj 
IFdPERT.GT.!l.A!IO.ISLIO.r.(Kl)CO  TO  90 
IFdPr.RT.LE.tOGO  TO  90 
C 

C  HOW  'C'  IS  THE  PCRTURRCO  MATRIX 

C  i:iv  A-  IHV  C+  FAC»d!IV  C)*(OA)*dllV  O* (FAC) •*2-( (MIV  C)» 

C  {nA))»»7.*dNV  C)+... 

C  HAPPRX-1  FIRST  2  TERMS  OF  THE  SERIES  FOR  IIIV  A  ARE  USED 

C  "2  FIRST  5  TERMS  OF  THE  SERIES  FOR  I’lV  A  ARE  USED 

C 

NAPPRX-2 

CALL  OPERT2(!l,MAX,!IAPPRX,FAC,C,DA,B) 

WRITE(G,llil) 

CALL  PI!MAT(II,MAX,B) 

GO  TO  2295 

C  APPLICATION  OF  PERTURBATI  0.1  METHOD  OVER - 

90  COIJTIilUE 

CALL  GKRPCT(!l,MAX,C,l!,nA) 

IF(PDIF-a0IF)CGG,65G,CG7 
GGG  CALL  C0RCT1(C,B,!I,MAX,5,1) 

GO  TO  GG8 

CG7  CALL  CORCT2(C,B,tl,MAX,5,l) 

668  CONTINUE 

IFdPERT.EO..O)GO  TO  2295 
IFdPERT.C’T.M.AHO.lSLID.EQ.DGO  TO  2291 
CALL  nPERTl(A,B,:;,MAX,FAC,  IPERT) 

2291  CONTINUE 

IFdPERT.LE.IOGO  TO  2295 
FFAC-FAC-l.On-OI) 

DO  91  K-1,11 
HRITC(G,300)K 

300  FORMAT(//,10X, 'VALUE  OF  K  -',12,/) 

DO  89  l=l,N 
DO  89  J-1,11 

89  Gd,J)-C(  l,J)-FFAC«DAd,J) 

|F(PDIF-QDIF)777,777,773 

777  CALL  CORCT1(G,B,N,MAX,5,0) 

GO  TO  779 

778  CALL  CORCT2(G,B,N,MAX,5',0) 

779  CONTINUE 
IF(K.LE.I»)FFAC=FFAC-10.0D0O 
IF{K.GT.'i.AND.K.LE.7  )FFAC»FFAC«1.773279D00 
I FCK.GE. 8 )FFAC"FFAC«1. 15478200 


6A 


i  B 


91 

2295 


580 


7S1 

782 

783 


92 


93 


94 


95 


80 

100 

101 

102 

103 

104 

105 
lOG 

107 

108 

109 

127 

150 

141 

140 


240 

249 

250 
550 
GOO 


COHTI.'UII' 

COliTItUli; 

ir(ii’!;nT. r.o..o., Min. ico:;cT.nQ. 0)00  to  733 
»ll!ITn(G,580) 

FORIlAT(5X,'ri:iAl.  lil!)  FOR  IJ'.l'UOVFr.CNT' ) 

F  RAC  •-!  .01)0 
II  IT- 11 

IF(PniF-QI)IF)7Sl,  731,782 
CALL  CORCTl(A,li,ll,llAX,IIIT,0) 

GO  TO  785 

CALL  C()RCT2(A,!i,ll,l!AX,l!n\,0) 

COllTll.'l)!; 

WI!ITL(0,109) 

CALL  PRMAT(H,MAX,i;) 

IHVERSL  AA  ■=  (tllVCRSC  Q)»(IMVCRGC  A)*(|:iVCRSn  P) 

CALL  I!F.0.UAT(!1.1*.AX,!5,C) 

IF(  ISP. CC.O.MIO.  ISO. Ffi. 0)00  TO  92 
CALL  llC(!UAT(!!,i;AX,n,l>\) 

IF{ iso.ro. DcALL  Dim .T(:j,:i,ti,i;\x,?\c,o,c,B,nA) 

I F  ( I p .  nr .  1 )  call  iif.i'.uat ( ii.iux, da, o 

ir(  ISP.r.O.  DCALL  lini!LT(:i,II,ll,MA;:,FAC,3,DA,P,C) 

CONTI  HUE 

IF(IPI!|:1T.C0.1)'IRITE(G,107) 

ir(  IPRINT. HO. DCALL  PRIlAT(:i,!UX,C) 

ir<PDir-0l)IF)95',93,54 

CONTMlur. 

CALL  i)Iiult(!I,:i,:i,i;ax,fac,o,a,b,g) 

IF(ISP.r.0.0)OO  TO  95 

CALL  niUILT(U,!l,i;,liAX,r\C,0,O,P,B) 

CALL  ni!iy(N,llAX,P) 

CAL L  D! HILT  ( II ,  N ,  11, MAX ,  F AC ,  0 ,  P,  B ,  G ) 

GO  TO  95 
coin  I  HUE 

CALL  Di;ULT(:i,l!,l!,MAX,FAC,0,B,A,G) 

1F(ISQ.!:Q.O)GO  to  95 

CA L  L  !)!  HU. T  ( 11 ,  :i ,  :i ,  1 1  AX ,  F A C ,  0 , 0,  G  ,  B ) 

CALL  DIIIV(!I,I;AX,0) 

CALL  PliULT(ll,:i,:!,!'AX,FAC,0,B,Q,G) 

CONTI  HUE 

IF<IPI{INT.r.O.D'.!RITr(G,103) 

I  F(  I  PHI  NT.  CO.  DCALL  PH:'.AT(N,I1AX,G) 

C A L  L  1 1C  R RO  H  ( ,  n AX ,  0  ,  C R RO R ) 
l)RITr.(G,24y) 
iF(iPi:in.cQ.o)sTOP 
COIITINUL 

FORMATC////////////////////) 

FORMAT (5X, 'DIAGONAL  SCALE  MATRIX,  P') 
rORMAT(5X,  '  IIJVCRSC  P  MAIRIX') 

F0R|;AT(5X, 'DIAGONAL  SCALE  IIATRIX,  O') 

FORMATC SX, ' INVCRSE  0  MATRIX') 

FORMAT(5X, 'MATRIX  A') 

FORIIATCSX,' IlIVCP.SC  A  MATRIX') 

rORMAT(5X,  '  IIIVERSC  OF  ORIGINAL  IIATRIX,  AA' ) 

FORIIATCSX, 'PRODUCT  OF  ORIGINAL  IIATRIX,  .AA,  AND  ITS  COMPUTED  IIIVERS 
IC) 

FORMATC /, lOX, ' IMPROVED  INVERSE  MATRIX') 

FORMATC5X, ' CROW)  SCALED  MATRIX  A') 

FORMATC 5 X, 'C«A+  r.PS«n  MATRIX') 

FORIIATCSX, '  INVCA)  "  INVCC)  ♦  CPS*(  )  ♦  CPS*»7  «C  )') 

FORMATC  lOX, 'MATRIX  DIMENSION  •' ,  I  2,//,  lOX,  '  ISP  ■= ' ,  I  2,/ ,1 9X,  '  I  SO  ■=' 
I,t2,/,10X,'IPE,RT  »',I2,/,10X,' ISLID  ,  I  2,/,  lOX,  '  ICOND  =',12, 
2/,10X, ' ICORCT  =',I2,//) 

FORMATC9I2) 

FORIIATCSX, 

1*.**. ***•****.**••*. .«..<*•*<*«***•**•••*..*•) 

FORMAT (3D25. 13) 

FORMATC/////, 5X, 'ORIGINAL  MATRIX,  AA’ ) 

FORMATC //,5X,' FAC" ',020. 11,/) 

STOP 

END 
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SUDUOUTiriC  GI'vRDCT(N,t)AX,X,Y,CCC) 

RCAL*3  X(M.\X,l),Y{HAX,l),A,Q,C,n,E,n!;T,CCC{a\X,l) 

IllTEfiCR  tlUM{2,2a) 

COtitiOtl  /  PR/ 1  PR  I  NT 
COMMOri  /  PROn/  pn  I F  ,  QD I  F 
00  G  1-1, (J 
00  G  J=l,il 
6  YCJ.n-XCJJ) 

A»1.0D0 
DO  G5  1=1,11 
0=0.000 
L=1 
M=l 
C 

C  FIND  LARGEST  ENTRY  A{U,M)  III  LOIIER  DIAGONAL  SUIIMATRIX 

C 

DO  18  J=l,ll 
DO  13  K=l,:i 

IF(I)AI1S(Y(K,J)).LC,D)G0  TO  18 
B-DAiJS(Y(R,J)) 

L=K 

H=J 

IS  CONTINUE 
C 

C  lUTERCIIANCin  ROWS 

C 

IFCL.LQ.DGO  TO  21) 

00  25  J  =  l,ll 
C=Y(L,U) 

Y(L,J)«Y(I,J) 

23  Y(I,J)=C 
C 

C  ItlTERClIAHGC  COLUMNS 

C 

24  IFdl.EQ.DGO  TO  29 
DO  28  J=1,N 
C=YtJ,M) 

Y{J,M)=Y(J,I) 

28  Y(J,I)=C 
C 

C  DF.GIN  SWEEP  COLUMNS  TO  THE  RIGHT 

C  ARRAYS  HUIUI,.)  ,:iU:i(2,.)  KEEP  RECORD 

C  OF  ROW  AND  COLUMN  INTERCHANGES 

C 

29  NUM{1,I)=L 
I1UM(2,I  )=M 
3=Y(I,I ) 

Y(I,I)=A 

DO  42  0=1, N 
IF(J.EQ.I)GO  TO  42 
C“-Y(I,J) 

Y(I,J)=O.ODO 
DO  41  K=1,N 
D=Y(i;,  I  )‘C 
E“Y(K,J)*t!+D 

lF(DAI'.S(r;).LT.1.0D-10»DABS;0))C«0,ODO 


C 


m  1 


i»  1 


Y{K,J)“E/A 

CONTINUE 

A-D 

RESTORE  COLUMNS 

00  50  I --2,11 

J=H+1-1 

K»NUH(2,J) 

IF<K.EQ.J)C.O  TO  52 
DO  51  L=a,ll 
C»Y(K,L) 

Y(K,L)*=Y(J,L) 

Y(J,L)=C 

K=llUIUl,d) 

RESTORE  ROWS 

IE(K.E(1.0)t!O  TO  53 
00  57  I -1,11 
C=Y{L,K) 

Y(L,K)“Y(L,0) 

Y<L,J)=C 

CONTINUE 

OET-'A 

WRITE(G,101)0ET 

F0I'.MAT(/1X, 'OET  OF  ORIG  MATRIX  IS  ',02lt.l7) 

OET-l.ODOO/OCT 

DO  111  1=1,11 

DO  111  J=1,N 

Y(I,J)=Y( l,J)*DET 

WRITE(G,102) 

FORIIATC/IX,  '  INVERSE  MATRIX  IS  •) 


CALL  PRHAT(II,I1AX,Y) 


001001=1, :i 
D0100J  =  1,!1 
CCC{|,J)=0,0000 

ooiooi;=i,:i 

IF(P0ir-Q0IF)Go,GG,G7 
CCC(I,J)=CCC(I,J)+X(I,K)»Y(K,J) 
GO  TO  100 

CCC(I,J)=CCC{I,J)+Y(I,X)*X(K.J) 

CONTINUE 


WRITE(G,103) 

I  F  (  I  PI!  1  NT .  EQ.  1  )CALL  PRMAT(!I,MAX,CCC) 

CALL  MERROR(ll,MAX,CCC,CR) 

F0RI1AT(/1X,' PRODUCT  OF  INVERSE  AND  ORIG  lUTRICES  '> 

RETURN 

END 
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SUBROUTINE  CORCTHX,  Y, N.MAXX  JIITCR,  I  PR) 

DlllCNSIO;i  X(r.AXX,l),Y(:!AXX,l),F{20,20),C(20,20),Z(20,23) 
COlUlON  /r)Ot./FRACI 
DOUBLE  PRECISION  X,Y,F,£,2 
DOUBLE  PRECISION  FRAC,ER1,ER2,CR 
DOUBLE  PRECISION  FRACI.FRSAV 
1/RITE(G,317)FI:ACI 
DO  200  ITER»1, NITER 
ER2“1.0ni0 
CR=ER2 
FRAC»FRACI 
EPa=O.ODOO 
C 

C  CALCULATE  ERROR  lUTRI X, F=X«Y- I 
C 

DO  210  1-1. N 
DO  210  J-l.N 
F(  l,J)=>0.0D00 
DO  200  K-l.N 

200  F('.J)nF(I.J)«XCI.K)*Y(K,J) 

IF(I .EQ.J)F( I,J)«F(I.J)-I.0n03 
ER1=ER1*F(I,.I)«F(I.J) 

210  CONTINUE 

CR1->DSQRT(ER1/(I1«II)) 

UI(ITC(G,320)ER1 

J20  FORItAT(//,10X,'ORIO|tlAL  RUSE  «  '.DIS.S,/) 

IF(  IPR.EQ.D'.IRITECG.SOO) 

IF(  IPR.LQ.DCALL  PRItATdl.n.F) 

C 

C  CALCULATE  IMPROVED  I NVCRSE, Y( IMPROVED)-Y-Y«F 
C 

DO  230  1  =  1, N 
DO  230  0  =  1,11 
E( l,0)=0.0D0 
DO  220  K=1,N 

220  E( l,J)  =  C(  l,0)*Y(  l,K)*F(K,0) 

250  CONTINUE 

IF(  IPR.EQ.l);lRITE(G,310) 

IF(  IPR.EQ.DCALL  PRUAT(N,12,C) 

255  CONTINUE 
ER2-EU 

CALL  MEOUAT(N,20.Y,Z) 

DO  2iiO  1  =  1, II 
DO  21(0  0  =  1, N 

21(0  Z(l,0)  =  Y(l,0)-FRAC*E(l,0) 

CALL  DMULT(N,N,':,20,FRAC,0,X,Z,F) 

CALL  MERROR(N,20,F,ER) 
ir(ER.OE.ER2)GO  TO  230 
FRSAV-FRAC 
FRAC-FRAC-O.snoo 
GO  TO  255 
230  CONTINUE 

DO  270  l-l,N 
DO  270  0=1, N 

270  Y(l,0)-Y(I,0)-FRSAV»E(I,0) 

290  CONTINUE 

300  FORHAT(5X,'F  MATRIX') 

310  FORMAT(5X,'E  MATRIX') 

317  FORMATCIOX,' INITIAL  FRAC-  '.Dll. 4) 

RETURN 

END 
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SUGROUT I  .'IE  COrvCT2  ( X ,  Y,  I!,  MAXX,  !l  11  ER,  t  PR) 

DlflE!lSIO:i  X(nAXX,l),Y(CAXX,l),r(;.),;0),EC>''^2P),2(20,?0) 
coiiiio:.'  /soL/ritACi 
DOUGLE  PRECISIOl;  X,Y,F,E,2 
DOUGLf  PKCCISIOIl  FRAv',E!a,CR2^ER 
OOUDLE  P.RECISIO:!  rR\CI,FRSAV 
WRITE(G,317)FRACI 
DO  290  1TEP.=>1,!!ITER 
CR2»1.0I)10 
ER-ER2 
FRAC-FRACl 
ER1>=0.0P00 
C 

C  CALCULATE  ERROR  MATP.I  X,  FoY*X- 1 
C 

no  210  1=1, :i 
no  210  0=1, M 
F(i,J)=o.onoo 

DO  200  K»l,:i 

200  F(l,J)»r(l,0)-*Y(l,K)*X(K,J) 

IF(I .C0.J)F( I,J)=F( l,J)-1.0n30 
CRi'EP.i+r(i,j)‘r(i,j) 

210  CORTUJUE 

ERl=nSQRT(ERl/(ll*«)) 

»RITE(G,320)CR1 

320  F0RMAT(//,10X, 'ORIGINAL  RMSE  ••  *,015.8,/) 

tF(IPR.Ea.l)WRITE(G,330) 

IFdPR.EO.DCALL  PR:;ATC),12,r) 

C 

C  CALCULATE  IHPROVCO  INVERSE, Y<  ll!PROVCO)=Y-F*Y 
C 

DO  230  l=l,:i 
DO  230  0--l,N 
E(  l,J)=0.0P0 
DO  223  K=1,N 

220  E(l,J)=E(l,J)*r(l,K)*Y(K,0) 

250  CONTIliUr 

IFUPR.EQ,l)l.'RITE(5,3ia) 

IF(IPR.E0..1)CALL  PR!!ATCI,12,E) 

255  CONTINUE 
ER2"CR 

CALL  ^•XQUAT(:!,20,Y,Z) 
no  2li0  1=1,11 
DO  21(0  3  =  1, N 

2!(0  Z(  l,J)=Y{l,JJ-FRAC*E<l,3) 

CALL  n.".ULT(N,I!,N,23,FRAC,0,Z,X,r) 

CALL  l!£RROR{N,20,F,CR) 

IF{ER.GE,ER2)GO  TO  730 
FRSAV=FRAC 
FRAC=FRAC*0.5D00 
GO  TO  255 
230  CONTI  NOE 

DO  270  1=1, N 
DO  270  3=1, N 

270  Y(I,3)=Y(I,3)-FRSAV«E(I,3) 

290  CONTINUE 

500  FORi’.AT(5X,'F  liATRIX') 

510  FORHAT(5X,'E  P.ATRi:;') 

317  FORtUTdOX,*  INITIAL  FRAC-  ',011,l|) 

RETURN 

END 
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SUDUOUT I  Il£  DPERTK  X,  Y, H.flAXX,  FAC,  I  PERT) 
dii;e:!sio.‘j  ;<(:iaaX,1),y<.-.axx,i),u{2D),v(20),f{23,2D) 
00U3LE  ppxcisio::  x,Y,u,y,F 
DOUGLE  PRECISION  rAC,ALPHA,nR 

coii:;o:j  /pr/ipr 


tPR=2 

IPR=1 

WRITr(G,21) 

CALL  0.';ULT(!!, 23, FAC  ,0,X,Y,F) 
tF(IPR.GE.l)i;RITE(3,95) 

I F ( I  PR . G E . 1 ) CALL  PRMATC N, 2 3, F ) 

CALL  MERROR(;!,20,F,Hn) 
write(6,320)e:; 

FORIIAK//, 10a, 'ORIGINAL  RMSE  =  *,015.8,/) 
CALCULATE  IliPROVEP  INVERSE 


ER=FAC 

FAC=FAC*):  ( I  PERT,  I  PERT) 

ALPUA=Y(IP£RT,II'ERT) 

ALP:iA=l .  300-rAC«ALPHA 
ALPi!A=rAC/ALPMA 

if(ipr.g:.i):ihitc(o,201)fac,alp«a 

fac=i;r 

00  150  1=1, N 
U(l)=Y(l,IPEr;T) 
v(I)=y(ipe:;t,  I) 
continue 

IF(IPn.CfJ.2);iRIT£(G,103)(U(l),l=l,N) 

IF(  IPR.r.S.2);.'RITr.(G,10S) 

I F  ( I  PR .  nn .  2  )'..'R  I TE  ( 0, 133  )  ( VC  I  ) ,  I =l,:i) 

DO  2^0  l=l,N 
no  21)0  J  =  1,N 
F(I,J)=U(  I)»AL?1.'A*VCJ) 

IF(  IPR.C5.1).IP.IT£(G,300) 

I F (  I  PR .  EQ .  1 )  CALL  PR! lAT ( N,  2 0, F  ) 

DO  250  1=1,11 
DO  253  J=1,N 
Y(1,J)=YC!,.I)*F(I,J) 

CALL  O.’lliLTCN, 20, FAC  .0,X,Y,F) 
IF(IPR.GE.l)w:;iTt(0.37) 

IF(  IPR.GE.DCALL  Pri:AT(;4,23,F) 

CALL  i;ekro!:(;i,20,f,er) 

WRITE (0,1 33) 

CALL  PR.':AT(II,|:AXX,Y) 

CONTINUE 

WRITE(6,91) 

^ FOR! ;AT ( 3 X,  '  APE RT :  A*  I  »V( C )  ' , / ) 
rOi::iAT(5X, 'APERT:  A*i:'.PKOVED  IMVCC)',/) 
FOi;:!AT(2:<,'i{i'25.1S,3X)) 

FORI-UC//) 

FOR;'.AT(/,10X,' l.'iPROVEn  INVERSE  MATRIX') 
FURIUTCSX,  'A  *EPS=',I)15.3,*ALPIIA=',IJ15.!),/) 
F0!’.i:AT(5X, 'F  MATRIX') 

RETURN 

END 
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SL'BROOTINR  DPERT2(N,ltAX,«APPRX,FAC,C,OA,0) 

DOUBLE  PRECISION  C(MAX,l),DA<tt\X,l),e{MAX,l),C(20,20) 
double  PRECtSIO.'l  FAC.POIF.QDIF 
COlUlOIJ  /PROD/POIF.QDIF 
CALL  GKRDCT(:J,tlAX,C,B,G) 

IF(P0IF-(J0IF)51,S1,S2 

51  CALL  CORCT1(C,B,IM!AX,N,0) 

GO  TO  53 

52  CALL  CORCT2(C,B,!l,:;AX,tl,0) 

53  CONTINUE 

CALL  nMULT(N,N,:!,!lAX,FAC,l,3,DA,G) 

CALL  DmiLT(:i/!,:!,MAX,FAC,3,r.,B,C) 

WRITECC.llS) 

CALL  PRCATdl.tlAX.C) 

IFCJAPPRX.Ea.DGO  TO  37 

CALL  D;iuLT(:i,iJ,:i,i:AX,FAC>o,r„G,nA) 

CALL  DMULTCi,. MAX, FAC,  0,P-A, C,C ) 

WRIT£<G,121) 

CALL  PRMAT(N,t!AX,G) 

87  COilTriUE 

DO  30  1-1,:.' 

DO  33  J-1,:: 

IF(IIArPi:X.E0..1)5(l,J)=B(l,J)+C(l,J) 

IF{:iAP?RX.EQ.2)5(l,J)»3(I.J)*C{l,J)«G(l,J) 

90  co:iti:jue 

118  F0PJ1AT(5X,'EPS«  MJV(C)  •D*IUV{C)‘) 

121  for:iat<5x,'eps*iiiv(C)  •')«i:;v(c)  ♦  cps*«2*( l'l•/(C)•D)••2 

RETURN 
END 


SUBROUTUiC  nTRA'JS(N,(WX,X) 

DOUBLE  PRECISIO;:  X(MAX,1),Y{20,20) 
C 

DO  1  i=i,;j 

DO  1  j=i.:j 

1  Y(I,J)=X(.I,I) 

DO  2  1=1,1! 

DO  2  J=1,N 

2  X(I,J)=Y(I,J) 

rctur:j 

END 


SUBfsOUT  I  tl  E  IIEO.U.AT  ( f '-'.X,  X,  Y) 

DOUBLE  PRECISION  X(I!AX,1),Y(."JVX,1) 

C 

C  EQUATE  MATRIX  Y  TO  MATRIX  X 

C 

DO  1  l=l,N 
DO  1  J-1,H 
1  Y(I,.I)»X(I,J) 

RETURN 

END 
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SURilOUTItin  l'SCALr<ll,f.AX,Tnacs,IS,»JF,X,Y) 
nOUIlLE  PRECISIO.I  X(I;AX,1),Y(!-.\X,1),R(20,:3) 

C0UR1.E  rtiUCISlO'J  T!iR£S,2LRO,3TOT,T'ir.,BT;:R,Y:-.\X,Y:'l!I,I'IF 

ROW  OR  COLUIRI  SCALiiir,  OF  i'ATRIX  X 
IS-1,  Ro;.'  GCALiriG 
IS=0,  COLU.-’.'l  SCAtllJG 

ZERO=1.0P-;o 

YtLAX=-1.0n*10 

Y!lli<=1.0R*10 

IF{  IS.EQ.  3)CAl.L  I1TRAHS(:J,,'1AX,X) 

00  30  1=1,-! 

T'IR'-SO.O 
00  10  J=l.!l 

tF(!)AnsCX(  !,J)).r.T.2ER0)n0  TO  5 

B(l,J)=-5«.3 

GO  TO  10 

co:!Ti;!tir 

t;(I,J)'DLO'.10(i}ADS(X(l,0))) 

IF(!i(l,J).GT,Ti:R)TilR=3(l,.0 

CORTHiOE 

ct:ir»tiir-tmrcs 

BTOT^O.0003 

ICT=0 

00  20  o»i>:: 

Y(I,J)=0.0030 
IF(C(I,J).LT.3T!!R)G0  TO  20 
ICTalCT+l 
CT0T=3T0T*3tl,J) 

C0«TI.'.'Un 

eTQT=>3TOT/lCT 

Y(l,l)-10.*«r.TOT 

1 F  ( l!TOT .  GT .  Yi-AX )  Y!  lAX  =»  RTOT 

I F  ( BTCT,  LT .  y;;  I  ;i )  Yi:  I  :J“8tot 

DO  50  j>=i.:i 

X(I,J)=X(I,J)/Y(I,I) 

COUTI.-IUE 

n!F=Yr.AX-yi;ir! 

I F  {  I S  .  Ed .  0 )  C AL  L  MTRA tis  C  H  .  f lAX  ,  X  ) 

RETUmi 

END 


SU3ROUTINC  r*!OLT(l'/I,L,IUX,S,IS,A,r.,C) 

DOUBLE  PRCCISIO.'J  A(;!AX,1 ), B(mX,l),C{mX,l  ),S,TCMP 

DOUBLE  PRECISION  f-UTRIX  liULTIPLIC.UIO;),  CsA*5«S 

DO  10  1=1, M 
00  10  J«1,L 
C(I,J)»0.0030 
DO  10  K=l,;i 

c(i,J)=S(i,J)+A(i,ro*n(K,j) 
IF(IS.EQ,1.A:JD.K.CQ,N)C(I,J)=S«C(I.J) 

10  CONTIHUE 
RETUR.'J 
EiiD 


72 


SpfSpliiM^ 


su[iROUTi:ir.  p'?;'\t(::.:'ax,\) 

Foi'.;;AT(2x,!i(n25.aG,5x)) 

FORf.AT(//) 

REAL*.*.  A((:\X,1) 

WRn£(5.105) 

K»((IJ-l)/'0  +  l 
DO  10  1=1, '! 

KRITE(5,107) 
no  15  1=1, t! 

J1  =  (L-1)*!»*1 

J2=J1*3 

iF(L.ra.r.)J2=-j 

WRITC(r.,103){A(l,J),J=Jl,J2) 

co:iTi'!i:c 
co:jti:)i;f 
WRITE  (0,1  oil) 
fori  !A1  (/////) 

FORWAK/) 

RETURI! 

END 


SUSROUTI iin  t<.£>:UOR(:!,:5AX,C,ERnOR) 
double  pr.ECisiO':  c(rAX,i),TDt,ATE:: 
DOUDLE  PREOISION  ERROR, ICR 
COiyiON  /IlC/AER 


ERROI'.=  0, 01)00 
AER«0,0n0 
DO  1  1=1, N 
DO  1  J=i,:i 
TEM-0(l,J) 

IF(I.EO.,J)Tr.M=C(l,O)-1.0Dne 

ATEK»DACS{TEM) 

I F  (  ATEi; .  GT ,  AE  P )  AER=A7Ell 
ERROR=Er.ROR*TE:;«TE:; 

CONTINUE 

ERROR-ERROR/ (:i*:i) 

ER30i;=nsar«T(  ERROR) 

WRITE(G,103)ERr,0R,AER 

FORILU(12X,’(."AT-l):  RHSE=' 

RETURN 

END 


,D11.I|,'  !1AX»',D11.!;) 


SUBROUTINE  OINV(N,IIAX,X) 

DOUBLE  PRECISION  X(|1AX,1) 

INVERSION  OF  DIAGONAL  IlATRIX,  X 

DO  1  l»l,!l 
DO  1  0=1, N 

IF{I.EQ.J)X(I,J)-I.0n00/X(l,0) 

RETURN 

END 


-  V^^'L 
'^2  = 


(CA) 


.‘5  .)  and  (C2)  may  bo  solved  simultaneously ,  to  yield  ll\o  desired 

r  function,  as  follows.  From  (01), 


V^-g,^K0,+C3)s  Yj 


e  “0  s 
^1  J 


V  » 
A 


(05) 


Vi+g2+(02+C3)«  -C.j.< 


Substitution  of  oqviation  (05)  Into  (02)  yields  the  ro(iuired  expression,  i.e.. 


IK 


S/ 


2,  0,  .s 

''  X-.C^s+l  J 

1.  M 


Yj+g^' (02+0j)«  Yj 


B3-0,s 


Y,+g2+(V‘:3)«  -‘y’ 


t0(>) 


gj-iy 


Y^+g^+O ,s 


Simp) if i';il ion  of  the  expression  for  H(s)  requires  unwieldy  algebiaie  manipulation, 
r.ni  w  'l  thoreiore  not  be  presented.  However,  it  can  be  shown  that  for  '/.^  '-ofil 
H(*^)  assumes  '  ne  following  form: 


H(s) 


K  s-fsli'^) 


(s+Oj)  (s-lo2^  ^4^ 


(07) 
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For  the  particular  transistor  circuit  shown  in  Fig.  Clb,  the  following 
parameter  values  are  assumed: 

=  O.Ollif  ^3  “  5  pf  83  ”  83  ~  AOmU  Zj^  =  1K52  (real) 

C2  =  )0.0pf  =  O.Olpf  g,  =  ]my  g2  =  0.5m<3 

Substituting  these  values  in  Equation  (C6)  and  performing  the  required  simplifications 

yields 

„(,) - 8JlCI^_.s^(s-8WI0(ipbj -  (es) 

(s+. 033 do'’) )  (s+. 080 (lo'’)  )  (s+25 . 2  (10*’) )  (s+1205 . 1  (10°) ) 

A  Bode  (corner) plot  of  Equation  (C8)  is  given  in  Fig.  C2  .  Inspection  of 
this  frequency  response  shows  that  the  circuit  is  broad-band  (i.e.,  it  poles  are 
separated  by  several  decades).  In  order  that  the  network  transfer  function  be 
identified  reliably,  the  spectrum  of  the  input  signal  must  contain  sufficient 
energy  concentrated  in  the  vicinity  of  the  network  poles  and  zeros  (for  excitation 
and  manifestation  of  these  critical  frequencies).  However,  it  is  not  convenient 
to  synthesize  —  and  realize  in  the  laboratory  —  an  input  having  such  characteristics. 

Therefore,  the  network  function  may  (and  in  this  case  will)  be  broken  into 
constituent  functions,  each  of  which  will  be  valid  for  a  particular  frequency  range. 

Inspection  of  Equation  (C8)  reveals  that  an  adequate  low-frequency  description 
of  H(s)  is 

2 

H  (s)  =  — (^1-vQ.Ai  - —  (Low-freq.)  (C9) 

"  (sf.033f;0  ')(s+.080(10  )) 

which  is  valid  up  to  1  Mr/s.  That  is.  Equation  (C9)  will  closely 

approximate  H(s)  for  radian  frequencies  below  approximately  10^  rad/sec.  This 
observation  can  be  seen  clearly  from  the  Bode  plot  in  Fig.  C2. 

Through  similar  considerations,  expressions  describing  the  mid-high  and  high 
frequency  characteristics  of  H(s)  are  obtained; 
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m  I 


IStt 


Equation  (CIO)  will  bo  valid  in  the 
b  9 

frequency  range  from  approximately  10  to  10  rad/sec,  while  Eeiuation  (i'll’) 

q 

will  be  valid  for  frequencies  from  10  rad/sec  onward. 

The  identification  technique  may  now  be  used  to  determine  models 
for  the  network  behavior  by  considering  each  of  the  three  regions  separately. 
Improvement  in  the  methodology  and  reliability  of  identification  of  broad-band 
networks  (systems)  is  being  investigated  under  a  new  research  task.  I'or 
example,  pre-filtering  the  output  data  in  order  to  isolate  the  variv  .is  frequenev 
regions  is  now  being  pursued. 


i*) 


Midband  gain  ■*  26.2  dB 

Fig.  C2.  Magnitude  characteristic  (Bode  plot) 
of  a  wide  bond  amplifier. 
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Our  approximate  description  of  the  low-frequency  behavior  of  H(s)  is 
given  in  eq.  (C9).  In  order  that  a  reliable  model  for  this  region  be  obtained, 
via  the  program  IGRAM,  several  factors  must  be  considered.  First,  a  careful 
choice  of  the  input  must  be  made  in  order  to  excite  the  low  frequency  modes  of  the 
system.  We  need  to  isolate  these  modes  of  the  response,  and  will  therefore 
use  an  input  signal  whose  spectral  content  is  concentrated  in  the  low 
frequency  region.  A  satisfactory  choice  is  a  single  triangular  pulse  of 
duration  125hsec.  This  signal  will  supply  sufficient  energy  to  the  low  frequency 
modes  and  relatively  small  amounts  to  the  higher  frequency  modes. 

Next,  we  must  decide  upon  a  sampling  interval,  A.  A  useful  rule-of-thumb 
in  making  this  choice  is  to  sample  at  a  frequency  f^  at  least  ten  times  the 
highest  frequency  of  interest.  For  the  system  under  consideration,  the  highest 
frequency  of  interest  is  0.013(10  )Hz*.  Thus,  a  sampling  interval  A  =  l/f^  = 

0.25  hsec  should  be  quite  adequate.  Notice  that  while  we  are  sampling  at  an 
adequate  rate  for  the  low  frequency  modes,  we  are  undersampling  the  high 
frequency  modes.  That  is,  the  system  as  a  whole  is  broad-band  and  we  are 
sampling  at  a  rate  suitable  only  for  the  low  frequency  portion.  Therefore, 
frequency  aliasing  can  be  expected  to  occur.  The  effect  of  this  aliasing,  however, 
(of  the  high  frequency  modes)  appears  as  evenly  distributed  noise  of  relatively 
small  power  spectrum  density. 

An  important,  but  less  obvious,  consideration  is  the  total  duration  of  the 
test  record  used  in  modeling.  Whenever  possible,  a  record  long  enough  to  have 
a  few  time  constants,  say  one  to  four,  of  the  slowest  mode  must  be  used.  Using 
this  criterion,  a  1000  point  record  (MPl  =  1000)  for  the  network  under  consideration 
should  suffice. 

ii)  Mid  to  High  frequency  Transistion 

Our  approximate  description  of  the  mid  to  high  frequency  transition  behavior 
of  H(s)  is  given  by  eq.  (CIO).  Considerations  similar  to  those  made  in  the  last 
section  yield  the  following  choices.  Realizing  that  a  narrowband  signal  must  be 


*  It  is  unrealistic  to  expect  that  the  design  or  test  engineer  know  the  exact 
frequencies  of  interest.  However,  it  is  assumed  that  he  has  some  idea  of 
the  critical  frequencies  of  the  system. 
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used  —  so  as  to  excite  only  the  mid-high  frequency  mode  (s=-25.2(10^)) ,  an 
exponentially  decaying  sinusoid  was  chosen  as  the  input  signal.  The  center 
frequency  of  this  input  lies  in  the  frequency  range  of  interest.  The  sampling 
interval  was  chosen  to  be  0.01  psec,  and  a  500-point  record  was  used  for  modeling 

In  modeling  this  region,  the  option  IBIAS  =  1  was  used.  The  reason  for 
this  choice  is  as  follows.  Due  to  the  low-frequency  modes,  a  transient  response 
will  appear  in  the  system  output  in  addition  to  the  desired  mid-frequency  respons 
However,  over  the  short  duration  of  our  record  (5psec)  this  slowly  varying 
transient  will  appear  relatively  constant,  resembling  a  d.c.  bias.  Ihe  option 
IBIAS  =  1  allows  the  program  to  separate  this  "bias"  and  hence  calculate  a 
more  reliable  model  for  the  mid-high  transition  range. 

iii)  High  Frequency  Region 

The  approximate  higli  frequency  description  of  H(s)  is  given  in  eq.  (Cll). 

The  input  signal  used  for  network  excitation  must  be  narrowband  (for  I'reviously 
mentioned  reasons).  Thus,  a  slowly  decaying  sinusoid  witl>  center  freqiieiuv  in 
the  critical  region  was  chosen.  Five  hundred  points  of  input-output  signals, 
with  a  sampling  interval  A  =  0.00025lisec,  were  used  for  modeling. 

Once  the  results  for  each  of  the  frequency  regions  have  been  obtained, 
they  may  be  used  to  synthesize  the  overall  network  response.  This  ^an  bi>  a.iu 
by  correctly  combining  the  model  descriptions  of  the  various  frequency  regiiuis. 
Details  of  sucli  a  syntliesis  will  not  be  discussed. 
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APPENDIX  D 

s-Domain  to  z-Domain  Conversion 

Samp led  SlRnal 

UTien  sampled  at  uniformly  spaced  time  instants  kA,  an  analog  signal 
x(t)  yields  a  numerical  sequence  C  -  ~  x(kA).  To  this 

CO 

numerical  sequence  we  can  associate  a  continuous- time  signal  x*(r)  =  Z 

k=-w 

Xj^6(t-kA),  called  the  (ideally-)sampled  signal.  If  the  original  signal  is 
bandlimited  by  1/2A  i'.z,  then  x(t)  can  be  recovered  from  x*(t)  through  low-pass 
filtering,  and  the  sampling  process  may  be  regarded  as  a  one-to-one  mapping. 

We  define  the  Laplace  transform  of  the  sampled  signal  in  the  customary  way; 
this  gives 


.X*(s)  = 


,  -sA.k 
x^(e  ) 


Now,  since  the  z  transform  of  C  =  {xj^}  is 


X(z)  =  I  (1^2) 

k— — CO 

we  make  the  extremely  interesting  observation  that 

X*(s)  =  X(z)  (D3) 

sA 

Note:  It  should  be  borne  in  mind  that  the  substitution  z=e  into 
X(z)  yields  the  Laplace  transform  of  x*(t),  not  of  x(t). 

Under  the  condition  of  bandlimitedness  (by  1/2A  Hz)  this 
substitution  yields  a  transform  that  agrees  with  X(s)  in  a 
suitable  neighborhood  of  s=C  in  the  s-plane. 

We  now  focus  attention  on  the  matter  of  conversion  of  transfer  functions 

from  s-domain  to  z-domain  and  vice-versa.  An  exhaustive  treatment  is  given 

in  reference  117].  Here,  we  summarize  three  of  the  most  widely  used  conversion 

techniques. 

1.  Logarithmic  Pole-Zero  Conversion 

sA  1 

This  technique  uses  the  relation  z  -  e  ,  or  s  =  ^  Lnz,  upon  the  poles 
and  zeros  of  the  function  under  consideration.  Thus 


s  TI  (s+b.) 
_ 

n 

Ti  (s+a  ) 
i=l  ^ 


iz-ir  (z-S  ) 

.  (n-£-m)  _ 1^1 _ 


Tt  (z-a  ) 
i=l  ^ 


i 

{ 


_ „.  "T 

E 


where 


-a^A 


a.  =  e 
1 


D5a 


6j  =  e 


-b.A 
1 


D5b 


2 .  Pu Ise-Invarlant  Conversion 

This  technique  has  the  merit  that  the  response  of  H(s)  to  an  input 


(t)  =  L  Xj_  p(t-kA),  where  p(t)  =  square  pulse  over  (0,A),  coincides 
k=-<» 


"  pulse 

with  the  response  of  K(2)  to  the  sequence  ^  many  cases  of  practical 

interest  i®  3n  excellent  approximation  to  x(t);  in  such  cases  this 

technique  of  conversion  promises  close  agreement  of  the  response  of  H{s)  to 
x(t)  and  of 
described  by 

n  b .  n  b  (1  -  e 

y  _ i _  u  / ,  ^  -  V  _  J _  _ 

-a.  A 


x(t)  and  of  H(z)  to  {x, ],  at  the  sampling  instants.  The  conversion  is 

K 


H(s)  =  I 


i=l  ®  ^  ^ 


HjCz) 


n 

r 

i=l 


-a  .A 

b^  (1  -  e  "  ) 


(D6) 


a.  (1  -  e  ^  z~^) 


3.  Impulse- Invariant  Conversion 

When  this  tecb.nique  is  used  the  response  of  H  (s)  to  x*(l)  coincides 
with  that  of  i'.vz)  to  =  (at  the  sampling  instants).  The  conversion  is 
described  by 


n 


b. 

1 


•'<«>  =  "  =  4.  . 

.  ,  s  +  a . 
1=1  1 


n  Ab . 
:^(z)  =  -  * 

i  =  l 


1  -  e 


-a  .A 
^  z-J 


(1)7) 


^ample:  Sampling  Inlerval  .A  =  5lis 


'iC 


-'WvV - 

Rj=10k 


—  AVW^ 

R2=100k 


H(s)  = 


+ 


;  X  10  _  _  . .  ^ 

s“  +  (1.2  X  10^)  s  +  (1  X  10') 


Cj  =  C,  =  lOOOpf 


0.025 


(s+9009.8)(s+110,990.2)  •'l "(z-0. 95595)  (z-0. 57410) 
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